-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.tex
1140 lines (1109 loc) · 71.4 KB
/
main.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[graybox]{svmult}
\usepackage{acronym}
\usepackage{amsmath,amssymb}
\usepackage[aboveskip=1pt,labelfont=bf,labelsep=period,justification=raggedright,singlelinecheck=off, font=small]{caption}
\usepackage{changepage}
\usepackage{cite}
\usepackage{nameref,hyperref}
\usepackage[utf8x]{inputenc}
\usepackage{newunicodechar}
\usepackage[right]{lineno}
\usepackage[nopatch=eqnum]{microtype}
\usepackage{xcolor}
\usepackage{multirow}
\usepackage{ctable} % for \specialrule command
\usepackage{multirow} % to combine cells vertically in a table
\usepackage{subcaption}
% ############ From Springer Template ###############
\usepackage{type1cm} % activate if the above 3 fonts are
% not available on your system
%
\usepackage{makeidx} % allows index generation
\usepackage{graphicx} % standard LaTeX graphics tool
% when including figure files
\usepackage{multicol} % used for the two-column index
\usepackage[bottom]{footmisc}% places footnotes at page bottom
\usepackage{newtxtext}
\usepackage[varvw]{newtxmath} % selects Times Roman as basic font
% ###################################################
\captionsetup[subfigure]{justification=centering}
\usepackage[outputdir=.]{minted}
\setminted{
linenos=true,
frame=single,
framerule=1pt,
obeytabs=true,
tabsize=4,
numbersep=-10pt
}
% Used to externally generate previews of the listings
\usepackage[tightpage]{preview}
\PreviewEnvironment{minted}
\bibliographystyle{unsrt}
\newcommand{\etal}{{\textit{et al. }}}
\newcommand{\mbx}{\mathbf{x}}
\newcommand{\mbu}{\mathbf{u}}
\newcommand{\mbp}{\mathbf{p}}
\newcommand{\mby}{\mathbf{y}}
\newcommand{\mbd}{\mathbf{d}}
\newcommand{\mbg}{\mathbf{g}}
\newcommand{\mbf}{\mathbf{f}}
\newcommand{\cmt}[1]{\textcolor{red} {[#1]}}
\providecommand{\keywords}[1]{\textbf{Keywords } #1}
% Create a new environment for code samples.
\newfloat{code}{tbp}{lop}
\floatname{code}{Code Sample}
% Create acronyms to use in text
\newacro{ode}[ODE]{Ordinary Differential Equation}
\newacro{oed}[OED]{optimal experimental design}
\newacro{fim}[FIM]{Fisher information matrix}
% Define unicode characters needed for output of simulation
\makeatletter
\verbatim@font
\newsavebox\v@sp
\sbox{\v@sp}{\ }
\newlength{\v@spwd}
\setlength{\v@spwd}{\wd\v@sp}
\newunicodechar{└}{\mbox{\kern0.5\v@spwd\vrule height 2ex depth -1ex width 0.2ex\kern-0.2ex\rule[1ex]{0.5\v@spwd}{0.2ex}}}
\newunicodechar{─}{\rule[1ex]{1\v@spwd}{0.2ex}}
\newunicodechar{┬}{\rule[1ex]{0.5\v@spwd}{0.2ex}\vrule height 1ex
depth 0.5ex width0.2ex\kern-0.2ex\rule[1ex]{0.5\v@spwd}{0.2ex}}
\newunicodechar{├}{\mbox{\kern0.5\v@spwd\vrule height 2ex depth 0.5ex width 0.2ex\kern-0.2ex\rule[1ex]{0.5\v@spwd}{0.2ex}}}
\newunicodechar{│}{\mbox{\kern0.5\v@spwd\vrule height 2ex depth 0.5ex width
0.2ex\kern-0.2ex \kern0.5\v@spwd}}
\makeatother
\normalfont
\begin{document}
% ########################################################################
% ########################################################################
\title*{Experimental design for predictive models in microbiology depending on environmental variables}
\titlerunning{Experimental Design in Predictive Microbiology}
\author{Polina Gaindrik, Jonas Pleyer, Daniel Heger and Christian Fleck}
\authorrunning{Polina Gaindrik \etal}
\institute{
Polina Gaindrik \at
\href{https://www.fdm.uni-freiburg.de/spatsysbio}{University of Freiburg, Freiburg Center for Data Analysis and Modeling},
Ernst-Zermelo-Str. 1, 79104 Freiburg, Germany,
\and
Jonas Pleyer \at
\href{https://www.fdm.uni-freiburg.de/spatsysbio}{University of Freiburg, Freiburg Center for Data Analysis and Modeling},
Ernst-Zermelo-Str. 1, 79104 Freiburg, Germany,
\and
Daniel Heger \at
\href{https://tsenso.com/en/}{tsenso GmbH, Johannesstr. 19, 70176 Stuttgart, Germany}
\and
Christian Fleck \at
\href{https://www.fdm.uni-freiburg.de/spatsysbio}{University of Freiburg, Freiburg Center for Data Analysis and Modeling},
Ernst-Zermelo-Str. 1, 79104 Freiburg, Germany,
\email{christian.fleck@fdm.uni-freiburg.de}
}
\maketitle
% ########################################################################
% ########################################################################
\abstract{
The aim of predictive microbiology is the provision of tools and methods for predicting the growth, survival, and death of microorganisms in different food matrices under a range of environmental conditions.
The parametrized mathematical models need to be calibrated using dedicated experimental data.
In order to efficiently plan experiments, model-based experimental design is used.
In this chapter, we explain model-based experimental design and provide step-by-step instructions for finding the optimal design using the well-known Baranyi-Roberts growth model as an example.
We provide the Python software {\it eDPM} for \ac{ode} -based models, such that the reader can apply model-based experimental design in their research context.
}
\keywords{Optimal experimental design, Parameter estimation, \ac{fim}, Identifiability, Uncertainty}
%
%
%
% ########################################################################
% ########################################################################
\section{Introduction}
Mathematical modeling is widely used to describe, understand and predict the behavior of living systems.
In particular, in the field of predictive microbiology, one can find a large variety of works that dwell on building models of different levels of complexity controlled by model parameters, e.g., to describe bacterial growth~\cite{bernaertsConceptsToolsPredictive2004}.
Predictive microbiology is a subfield of food microbiology based on the premise that the behavior of microorganisms can be described by the use of mathematical models accounting for the effects of various factors such as temperature, pH, water activity, and the presence of preservatives, on the growth and survival of microorganisms.
These models can then be used to predict the behavior of microorganisms under different conditions and to evaluate the effectiveness of various control measures, such as refrigeration, heat treatment, or the addition of preservatives.
Predictive microbiology has many applications in food safety, as it can be used to assess the risk of foodborne illness, to design safe food processing and storage practices, and to develop new food products with extended shelf-life.
For successful predictions it is essential that the model parameters can be estimated from experimental data.
However, due to measurement noise the estimates are accompanied by some uncertainties.
That gives rise to a set of important questions:
Can all parameters of the given model structure be estimated?
What should be measured and when? What are the confidence intervals for the parameters?
How can the experimental effort be minimized?
Answering these questions equates to finding the \ac{oed}.
This results in defining optimized experimental conditions and/or measurement times which allows for a reduction of the experimental load without loss of information~\cite{derlindenImpactExperimentDesign2013, balsa-cantoe.bangaj.r.COMPUTINGOPTIMALDYNAMIC2008}.
In general, experimental design can be also used for model discrimination~\cite{kreutzSystemsBiology2009, stamatiOptimalExperimentalDesign2016}.
However, in this chapter we focus on the design procedure developed to increase parameter estimation precision.
Comprehensive reviews on can be found here:~\cite{atkinsonDevelopmentsDesignExperiments1982, franceschiniModelbasedDesignExperiments2008,sunParameterEstimation2011,vilasPredictiveFood2016}.\\
The whole parameter estimation process starts with choosing a model structure and defining observable quantities which the researcher is interested in (see Fig.~\ref{Fig1}).
Once this is chosen, a first parameter set is selected based on literature review or prior data.
In order to fine-tune the experimental design setup to the given use-case, we need to choose an objective function that will be optimized.
Depending on the goal, a researcher faces an important choice; which of the several sometimes contradicting objectives should be chosen for a particular case?
For example, is it more desirable to have precise knowledge in a chosen set of parameters and less information about the remaining ones?
What is the best balance between experimental effort and precision?
Using multi-objective approaches, some attempts were made to answer these questions and to improve the experimental design by combining several objectives~\cite{telenOptimalExperimentDesign2012, logistRobustMultiobjectiveOptimal2011}.
For models that depend on external cues like temperature, it is {\it a priori} not clear whether constant variable conditions are more efficient for parameter estimation~\cite{versyckIntroducingOptimal1999,garciaQualityShelflifePrediction2015}.
Using the initial guess of parameters a first \ac{oed} is determined taking into account specific constraints, such as maximum number of measurements, measurement times, etc.
Next, the designed experiments need to be performed and this data results into a new estimation of the parameters.
If the confidence intervals of the parameters are sufficiently small, the scheme ends, otherwise one uses the new estimates as the starting point for the next experimental design.
The process can be repeated several times to increase the precision of the parameter estimates until the desired accuracy is achieved.
In order to test the scheme one can also perform numerical experiments and obtain {\it in-silico} data.
%
%
%
% ########################################################################
% ########################################################################
\section{Materials}
There are several software packages available for model-based experimental design~\cite{balsa-canto_amigo2_2016, zhang_optimal_2018, busetto_near-optimal_2013}.
These packages comprise numerous tools and consequently require some time to understand how to use them.
We developed {\it eDPM} (Experimental Design for Predictive Microbiology), which focuses on ease of use.
It is a set of tools to calculate the sensitivity and Fisher information matrix and do parameter space exploration in order to find optimal results.
We require a working installation of the popular scripting language Python $\geq3.7$~\cite{rossumPythonLanguageReference2010}.
For installation instructions, please refer to the website \href{https://www.python.org/downloads/}{python.org}.
In addition, we expect users to be able to write, execute and display output of scripts.
This tutorial can also be followed using Jupyter Notebook~\cite{jupyterteamJupyterNotebook}.
%
Users can obtain it by installing {\it eDPM} from \href{https://pypi.org/project/edpm/0.0.1/}{pypi.org}.
The python website has guides for installing packages \href{https://packaging.python.org/en/latest/tutorials/installing-packages/}{packaging.python.org/}.
Most Unix-Based systems (e.g., GNU/Linux) and Mac-OS can use \mintinline[bgcolor=white,style=emacs]{bash}{pip}
\begin{minted}[bgcolor=white, linenos=false]{bash}
pip install eDPM
\end{minted}
or \mintinline[bgcolor=white,style=emacs]{bash}{conda} to install the desired package.
\begin{minted}[bgcolor=white, linenos=false]{bash}
conda install eDPM
\end{minted}
The \href{https://spatial-systems-biology-freiburg.github.io/eDPM/}{documentation} of the package is continuously updated.
After the installation of the package is complete, we open a file with a text editor of choice and simply start writing code.
We begin by writing a so-called she-bang which is responsible to signalize that this file is to be executed with python.
Afterwards, we import the needed packages.
\begin{code}[h]
\begin{minted}[firstnumber=1]{python}
#!/usr/bin/env python3
import numpy as np
from eDPM import *
\end{minted}
\caption{Import statements to use {\it eDPM}}
\label{code:import_statements}
\end{code}
This will serve as the starting point for our script.
In the following, we will append more code to it and utilize the methods developed in {\it eDPM}~\cite{edpm2023}.
The Code Samples containing line numbers correspond to one script that defines the system subjected to optimization.
The line numbers indicate the order of the code.
%
%
% ########################################################################
\section{Methods}
\subsection{Model Formulation}
\subsubsection{Theory}
As a starting point it is necessary to define the mathematical model.
We restrict our discussion to biological systems which can be described by a system of \aclp{ode}:
\begin{equation}
\begin{cases}
\dot \mbx (t) = \mbf(t, \mbx, \mbu, \mbp) \\
\mbx (t_0) = \mbx_0.
\label{eq:ode_general}
\end{cases}
\end{equation}
Here $\mbx = (x_1, x_2, ..., x_n)$ is a vector of state variables of the system with initial condition $\mbx_0$, $t$ is a time, $\mbu$ is a vector of externally controlled inputs into the system (e.g., supply of glucose or external temperature) and $\mbp$ are the parameters of the system.
We assume that a subset of the parameters $\mbp$ is unknown and should be estimated from data.
Predominantly, the state variables cannot be directly observed, but rather the observables are functions of the state variables:
\begin{equation}
\mby (t) = \mbg(t, \mbx (t), \mbu, \mbp) + \epsilon (t),
\label{eq:observ_general}
\end{equation}
where the function $\mbg$ is an observation or output function, and $\epsilon$ is a measurement noise.
The observational noise is often assumed to be an independently distributed Gaussian noise with zero mean and variance $\sigma^2$: $\epsilon (t) \sim N(0, \sigma^2)$, $\forall t$.
In this tutorial, we demonstrate the usage of experimental design on the widely employed mathematical model developed by Baranyi and Roberts~\cite{baranyiDynamicApproach1994}, which was devised to describe bacterial growth.
The model introduces two state variables $\mbx = (x_1, x_2)$, where $x_1(t)$ denotes the cell concentration of a bacterial population at the time $t$
and $x_2(t)$ is the quantity defining the proportion of the growth rate specified by the environment, e.g., a limiting nutrient critical for bacterial growth:
\begin{equation}
\begin{cases}
\dot x_1(t) = f_1(x_1, x_2) = \frac{x_2(t)}{x_2(t) + 1} \mu^\text{max} \big(1 - \frac{x_1(t)}{x_1^\text{max}}\big) x_1(t)\\
\dot x_2(t) = f_2(x_1, x_2) = \mu^\text{max} x_2(t).
\label{eq:ode_BaranyiRoberts}
\end{cases}
\end{equation}
Here $\mu^\text{max}$ determines the maximum growth rate, and $x_1^\text{max}$ is the maximal bacteria concentration due to environmental constraints.
The condition $x_2(0)$ allows to quantify the initial physiological state of the cells and, hence, the process of adjustment (lag-phase)~\cite{grijspeerdt_estimating_1999}.
To account for the influence of the temperature on the activity of the model, we will use the 'square root' or Ratkowsky-type model for the maximum growth rate~\cite{ratkowsky_relationship_1982}
\begin{equation}
\sqrt{\mu^\text{max}} = b (T - T^\text{min}),
\label{eq:RatkowskyModel}
\end{equation}
where $b$ is the regression coefficient, and $T^\text{min}$ is the minimum temperature at which the growth can occur.
As the observable we choose the bacterial count, i.e.:
\begin{equation}
y(t) = x_1(t)+\epsilon(t),
\label{eq:observable}
\end{equation}
a common choice in predictive microbiology.\\
The Eqs.~\ref{eq:ode_BaranyiRoberts} - \ref{eq:observable} fully define the system.
Here $x_1^\text{max}, b, T^\text{min}$ are the parameters that we estimate using observational data $y$ at measurement times $t$, and temperature $T$ is an input of the system.
Based on this model, we would like to optimize the choice of measurement times as well as temperatures (inputs) of the system to find the \acl{oed}.
%
% ########################################################################
\subsubsection{Code}
In order to be able to solve the equations numerically, we first need to define the \ac{ode} system described by Eq.~\ref{eq:ode_general} in python.
Next, we determine all numerical values present in the system.
We distinguish between time points $t_i$ at which the result of the \ac{ode} is evaluated, inputs $\mbu$, which alter the behavior of our system (for example temperature, humidity, etc.), parameters $\mbp$, which are the quantities that we want to estimate, and other arguments needed to solve the \ac{ode}.
Table \ref{Table1} provides an overview of these quantities and gives corresponding names in the code.
In the following, we will step-by-step explain how to specify all variables needed by using the example of the Baranyi-Roberts model.\\
%
First we define the right-hand side of the \acp{ode}, i.e., we need to implement the Baranyi-Roberts model (Eqs.~\ref{eq:ode_BaranyiRoberts} - \ref{eq:RatkowskyModel}) into a function as can be seen in Code Sample~\ref{code:code_baranyi_roberts_ode_fun}.
\begin{code}[h]
\begin{minted}[firstnumber=6]{python}
# The function name can be chosen freely, but the order
# of required function arguments is fixed by the definition.
def baranyi_roberts_ode(t, x, u, p, ode_args):
# Unpack the input vectors x,u,p for easy access
# of their components
(x1, x2) = x
(Temp, ) = u
(x_max, b, Temp_min) = p
# Calculate the maximum growth rate
mu_max = b**2 * (Temp - Temp_min)**2
# Calculate the right hand side of the ODEs, store
# it in a list [ ... ] and return it.
return [
mu_max * (x2/(x2 + 1)) * (1 - x1/x_max) * x1,
mu_max * x2
]
\end{minted}
\caption{
Definition of the Baranyi-Roberts \ac{ode} model.
}
\label{code:code_baranyi_roberts_ode_fun}
\end{code}
%
% ########################################################################
\subsection{Parameter Estimation}
After defining the model, the user needs to provide an initial parameter set.
This could be chosen from the literature or estimated from previously gathered experimental data.
It is common to assume that the observational or measurement noise is Gaussian white noise (e.g., no temporal correlations) with zero mean and variance $\sigma^2(t)$: $\epsilon(t) \sim N(0, \sigma^2(t))$, $\forall t$~\cite{kreutzProfileLikelihood2013}.
In this case the logarithm of the likelihood function (log-likelihood) is given by Eq.~\ref{eq:likelihood_Gaussian}:
\begin{equation}
\ln L(\mbp) \propto - \sum_{i}\frac{ \big(\mathbf{g}^{i}(\mbp) - \mbd^{i}\big)^2}{2 \sigma_{i}^2}.
\label{eq:likelihood_Gaussian}
\end{equation}
Here we introduced the following shorthand notations:
\begin{itemize}
\item $\mbd^{i}$: data measured at time $t=t_i$
\item $\mathbf{g}^{i}(\mbp)$: the observation function depending on the parameter vector $\mbp$ evaluated at time $t=t_i$
\item $\sigma_{i}^2$: variance of the observational noise at time $t=t_i$.
\end{itemize}
The method of maximum likelihood consists in searching for the parameter vector $\mbp$ which maximizes the (log-)likelihood~\cite{gaborRobustEfficient2015}.
%
%
% ########################################################################
\subsection{Experimental Design}
Following our definition of the model and setting the initial parameter vector $\mbp_0$ we can proceed with experimental design.
In essence, experimental design comprises maximizing an objective function depending on the model, the initial parameter vector $\mbp_0$, external input $\mbu$ into the system, and a set of discrete observation times $t_i$ at which the potential measurements should be performed.
The objective function needs to quantify the information the observation $\mby$ has about the parameter vector $\mbp$.
A common choice here is constructing objective functions based on the \acf{fim}~\cite{lyTutorialFisher2017}.
According to the so-called 'Cramer-Rao inequality', the Fisher information is inversely proportional to the minimal squared estimation error~\cite{friedenExploratoryData2010}.
This relation justifies using the FIM to improve our experimental design.
In the following we explain one of the ways to calculate the \ac{fim} for an \ac{ode} system (Eq.~\ref{eq:ode_general}) and the observable function (Eq.~\ref{eq:observ_general})~\cite{lyTutorialFisher2017}.
%
%
% ########################################################################
\subsubsection{Sensitivity Calculation}
An easy way to calculate the \ac{fim} uses local sensitivities of the observables $\mby$, that are in turn calculated from local sensitivities of the internal variables $\mbx$~\cite{versyckIntroducingOptimal1999, banks_generalized_2010}.
Assume that functions $\mbf$ and $\mbg$ are differentiable functions with respect to the state variables $\mbx$ and parameters $\mbp$.
The local sensitivities of $i$-component of the vector $\mbx$ w.r.t. $j$-component of the parameter vector are defined by $s^x_{ij} = (\mathrm{d} x_i / \mathrm{d} p_j)$.
These can be calculated using an enhanced system of the \acp{ode}:
\begin{equation}
\begin{cases}
\dot x_i (t) = f_i(t, \mbx, \mbu, \mbp)\\
\dot s^x_{ij} = \sum_k \frac{\partial f_i}{\partial x_k} s^x_{kj} + \frac{\partial f_i}{\partial p_j}.
\label{eq:ode_and_sensitiv}
\end{cases}
\end{equation}
For the \ac{fim} we need the sensitivities of observable functions $(\mathrm{d} y_i / \mathrm{d} p_j) = s_{ij}$.
We determine these at the certain times $t_m$ and input $u_n$ using the solutions of $s^x_{ij}$:
\begin{equation}
s_{ij} (t_m, u_n) = \sum_k \frac{\partial g_i}{\partial x_k}\bigg|_{t_m, u_n} s_{kj}^x (t_m, u_n) + \frac{\partial g_i}{\partial p_j}\bigg|_{t_m, u_n}.
\label{eq:observ_sensitivities}
\end{equation}
\\
These sensitivities are the elements of the sensitivity matrix~\cite{stigterObservabilityComplex2017}.
For example, in case of two observables $\mby = (y_1, y_2)$, two different inputs $\mbu = (u_1, u_2)$, $N$ different measurement times and $N_p$ parameters, sensitivity matrix reads:
%
\begin{equation}
S =
\begin{pmatrix}
s_{11} (t_1, u_1) & ... & s_{1 N_p}(t_1, u_1) \\
\vdots & & \vdots \\
s_{11} (t_{N}, u_1) & ... & s_{1 N_p} (t_{N}, u_1)\\
s_{11} (t_1, u_2) & ... & s_{1 N_p}(t_1, u_2) \\
\vdots & & \vdots \\
s_{11} (t_N, u_2) & ... & s_{1 N_p} (t_N, u_2)\\
s_{21} (t_1, u_1) & ... & s_{2 N_p}(t_1, u_1) \\
\vdots & & \vdots \\
s_{21} (t_{N}, u_1) & ... & s_{2 N_p} (t_{N}, u_1)\\
s_{21} (t_1, u_2) & ... & s_{2 N_p}(t_1, u_2) \\
\vdots & & \vdots \\
s_{21} (t_N, u_2) & ... & s_{2 N_p} (t_N, u_2)
\end{pmatrix}
\label{eq:sens_matrix}
\end{equation}
This matrix will be used to directly calculate the \ac{fim}:
\begin{equation}
F = S^T Q^{-1} S.
\end{equation}
Here, $Q$ is the covariance matrix of measurement error~\cite{cintron-arias_sensitivity_2009, balsa-cantoComputationalProcedures2008}.
If the measurements are independent, then only the diagonal elements of the matrix are non-zero:
\begin{equation}
Q =
\begin{pmatrix}
\sigma_{1}^2(t_1, u_1) & 0 & 0 & \dots & 0 \\
0 & \sigma_{1}^2(t_2, u_1) & 0 & \dots & 0 \\
0 & 0 & \sigma_{1}^2(t_1, u_2) & \dots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \dots & \sigma_{2}^2(t_N, u_2)
\end{pmatrix},
\label{eq:covar_matrix}
\end{equation}
where $\sigma_{i}^2 (t_m, u_n)$ is the error of the measurements of observable $y_i$ at time point $t_m$ with input $u_n$.
The error can be either estimated from data or approximated by some error model.
For example, one may consider that the error contribution consists of an absolute part that stays constant and a relative part that is proportional to the observable values.
Then the diagonal elements of the covariance matrix take form:
\begin{equation}
\label{eq:error_model}
\sigma_{i}^2 (t_m, u_n) = \gamma_\text{abs} + \gamma_\text{rel} \cdot y_i(t_m, u_n),
\end{equation}
where $\gamma_\text{abs}$ and $\gamma_\text{rel}$ are the coefficients determining the absolute and the relative error contribution, respectively.
\\
In the case the parameters are of very different scale (e.g., on the second and on the day scale), it may be advisable to use the relative or normalized sensitivities to improve not the absolute but the relative accuracy of the parameter estimates:
\begin{equation}
\tilde{s}_{ij} (t_m, u_n) =\frac{\mathrm{d}\ln(y_i)}{\mathrm{d}\ln(p_j)} = \frac{\mathrm{d} y_i}{\mathrm{d} p_j} \frac{p_j}{y_i}\bigg|_{t_m, u_n}.
\label{eq:relat_sensitivities}
\end{equation}
In this instance one needs to use relative measurements errors for the matrix $Q$:
\begin{equation}
\tilde{\sigma}_{i}^2 (t_m, u_n) = \frac{\sigma_{i}^2 (t_m, u_n)}{y^2_i(t_m, u_n)}.
\end{equation}
%
% ########################################################################
\subsubsection{Numerical Sensitivity Calculation}
For the calculation of the sensitivities (Eq.~\ref{eq:ode_and_sensitiv}) we need to supply the derivatives of $f$ with respect to the parameters and state variables.
We define $\mbx=(x_1,x_2,x_3)$ and $\mbp=(p_1,p_2,p_3)$ with $p_1=x_1^\text{max}$, $p_2=b$ and $p_3=T^\text{min}$.
Firstly, we calculate the mathematical derivatives $\partial f_i/\partial p_j$ and $\partial f_i/\partial x_j$ and then implement the corresponding functions.
The first component of the Baranyi-Roberts model reads:
\begin{equation}
\label{eq:rob_bar_model}
\dot x_1(t) = f_1(t, \mbx, \mbu, \mbp) = \mu^\text{max} \frac{x_2(t)}{x_2(t) + 1} \left(1 - \frac{x_1(t)}{x_1^\text{max}}\right) x_1(t),
\end{equation}
where $\sqrt{\mu^\text{max}}=b(T-T^\text{min})$.
It follows:
\begin{alignat}{6}
\frac{\partial\mu^\text{max}}{\partial p_1} &= 0\\
\frac{\partial\mu^\text{max}}{\partial p_2} &= 2b(T-T^\text{min})^2 &=& &&\frac{2\mu^\text{max}}{b}\\
\frac{\partial\mu^\text{max}}{\partial p_3} &= -2b^2(T-T^\text{min}) &=& &-&\frac{2\mu^\text{max}}{T-T^\text{min}}
\end{alignat}
and consequently:
\begin{alignat}{7}
\frac{\partial f_1}{\partial p_1}(t) &= \mu^\text{max} &&\frac{x_2(t)}{x_2(t) + 1} &\frac{x_1(t)}{\left(x_1^\text{max}\right)^2} &&x_1(t)\\
\frac{\partial f_1}{\partial p_2}(t) &= \frac{2\mu^\text{max}}{b} &&\frac{x_2(t)}{x_2(t) + 1} &\left(1 - \frac{x_1(t)}{x_1^\text{max}}\right) &&x_1(t)\\
\frac{\partial f_1}{\partial p_3}(t) &= -\frac{2\mu^\text{max}}{T-T^\text{min}} &&\frac{x_2(t)}{x_2(t) + 1} &\left(1 - \frac{x_1(t)}{x_1^\text{max}}\right) &&x_1(t).
\end{alignat}
Similarly, the derivatives $\partial f/\partial x_i$ are given by:
\begin{alignat}{7}
\frac{\partial f_1}{\partial x_1}(t) &= \mu^\text{max} &\frac{x_2(t)}{x_2(t) + 1} && \left(1-2\frac{x_1(t)}{x_1^\text{max}}\right) &&\\
\frac{\partial f_1}{\partial x_2}(t) &= \mu^\text{max} &\frac{1}{(x_2(t)+1)^2} && \left(1 - \frac{x_1(t)}{x_1^\text{max}}\right) &&x_1(t).
\end{alignat}
The readers are encouraged to calculate the components $\partial f_2/\partial p_i$ and $\partial f_2/\partial x_i$ as an exercise.
The resulting implemented functions in python can be seen in Code Sample~\ref{code:code_baranyi_roberts_ode_fun_derivatives}.
\begin{code}[H]
\begin{minted}[firstnumber=25]{python}
def ode_dfdp(t, x, u, p, ode_args):
(x1, x2) = x
(Temp, ) = u
(x_max, b, Temp_min) = p
mu_max = b**2 * (Temp - Temp_min)**2
return [
[
mu_max * (x2/(x2 + 1)) * (x1/x_max)**2,
2 * b * (Temp - Temp_min)**2 * (x2/(x2 + 1))
* (1 - x1/x_max)*x1,
-2 * b**2 * (Temp - Temp_min) * (x2/(x2 + 1))
* (1 - x1/x_max)*x1
],
[
0,
2 * b * (Temp - Temp_min)**2 * x2,
-2 * b**2 * (Temp - Temp_min) * x2
]
]
def ode_dfdx(t, x, u, p, ode_args):
(x1, x2) = x
(Temp, ) = u
(x_max, b, Temp_min) = p
mu_max = b**2 * (Temp - Temp_min)**2
return [
[
mu_max * (x2/(x2 + 1)) * (1 - 2*x1/x_max),
mu_max * 1/(x2 + 1)**2 * (1 - x1/x_max)*x1
],
[
0,
mu_max
]
]
\end{minted}
\caption{Derivatives of the function $f$ of the Baranyi-Roberts model \ac{ode}.}
\label{code:code_baranyi_roberts_ode_fun_derivatives}
\end{code}
%
% ########################################################################
\subsubsection{Define numerical values}
After the definition of the structure of the \acp{ode}, we need to specify numerical values.
It is good practice to gather such definitions in the \mintinline[bgcolor=white,style=emacs]{python}{__main__} method of the python program as it was done in Code Sample~\ref{code:code_write_main_function_start}.
We start by defining the parameters of the \acp{ode} (line 64) and afterwards the initial values (line 67).
Next, we constrain the optimization routine to find the best possible time points in the interval $t_i\in\left[t_\text{low},t_\text{high}\right]$.
For the purpose we write the times as a dictionary with entries \mintinline[bgcolor=white,style=emacs]{python}{{"lb":t_low, "ub":t_high, "n":n_times}} where \mintinline[bgcolor=white,style=emacs]{python}{"lb"} and \mintinline[bgcolor=white,style=emacs]{python}{"ub"} are the lower and the upper bound while \mintinline[bgcolor=white,style=emacs]{python}{n} describes the number of discrete time points at which we want to sample the system (line 70).
In contrast, if we wanted to specify explicit values for the sampling points, we would have needed to supply a list of time values or a \mintinline[bgcolor=white,style=emacs]{python}{np.ndarray} directly.
One can see this approach in lines 78-80 of the code example.
Here, we supply an \mintinline[bgcolor=white,style=emacs]{python}{np.ndarray} with explicit values, thus fixing the sampling points for the optimization routine.
\begin{code}[h]
\begin{minted}[firstnumber=62]{python}
if __name__ == "__main__":
# Define parameters
p = (2e4, 0.02, -5.5)
# Define initial conditions
x0 = np.array([50, 1])
# Define interval and number of sampling points for times
times = {"lb":0.0, "ub":100.0, "n":2}
# Define explicit temperature points
Temp_low = 2.0
Temp_high = 12.0
n_Temp = 2
# Summarize all input definitions in list (only temperatures)
inputs = [
np.linspace(Temp_low, Temp_high, n_Temp)
]
\end{minted}
\caption{
The main function will encompass every step in our experimental design approach.
In the beginning, we insert the actual values for our model definition.
}
\label{code:code_write_main_function_start}
\end{code}
\subsubsection{Defining Explicit Values and Sampling}
There are two distinct ways to define inputs of the model.
First, one can fix (multiple) explicit values for an input variable
\begin{minted}[linenos=false]{python}
inputs = [
np.linspace(Temp_low, Temp_high, n_Temp)
]
# >>> inputs
# [array([2., 12.])]
\end{minted}
so that experiment will then be solved for only these values without optimizing them.
Secondly, one can also optimize inputs of the system.
For example, suppose we want to optimize the temperature at which the experiments are performed, which means we let the optimization algorithm pick the optimal temperature points such that the information gathered from the system is maximized.
To sample in the interval \mintinline[bgcolor=white,style=emacs]{bash}{(Temp_low, Temp_high)} with \mintinline[bgcolor=white,style=emacs]{bash}{n_Temp} points the reader can simply write:
\begin{minted}[linenos=false]{python}
inputs = [
{"lb":Temp_low, "ub":Temp_high, "n":n_Temp}
]
# >>> inputs
# [{'lb': 2.0, 'ub': 12.0, 'n': 2}]
\end{minted}
The difference between choosing explicit values and specifying a sampling range may be subtle at first glance, but in turn allows to very easily switch between fixed values and optimized sampling points.\\
%
The different variables can be treated individually as needed.
Suppose, we have a system with temperatures and humidity as input variables.
We could decide to have the temperature optimized, but fix the humidity explicitly, because experimentally one can change the temperature continuously (or in discrete steps) but is restricted regarding the humidity.
In our code, we simply would mix explicit and sampled definitions for individual inputs:
\begin{minted}[linenos=false]{python}
inputs = [
# These values will be sampled
{"lb":Temp_low, "ub":Temp_high, "n":n_Temp},
# These are fixed
np.array([0.4, 0.45, 0.5, 0.55])
]
\end{minted}
%
% ########################################################################
\subsubsection{Define Model}
After we have decided on the numerical values for our model, we need to put everything together.
The \mintinline[bgcolor=white,style=emacs]{python}{FisherModel} class serves as the entry point.
Here, we simply supply every previously made definition of variables and methods (see \textbf{Note~\ref{note-order-arguments-python}}).
\begin{code}[H]
\begin{minted}[firstnumber=82]{python}
# Create the FisherModel which serves as the entry point
# for the solving and optimization algorithms
fsm = FisherModel(
ode_x0=x0,
ode_t0=0.0,
ode_fun=baranyi_roberts_ode,
ode_dfdx=ode_dfdx,
ode_dfdp=ode_dfdp,
times=times,
inputs=inputs,
parameters=p
)
\end{minted}
\caption{Define the full model.}
\label{code:model_definition}
\end{code}
\noindent For a full list of optional arguments, we refer to the \href{https://spatial-systems-biology-freiburg.github.io/eDPM/}{documentation} of the package.
%
% ########################################################################
%
% ########################################################################
\subsubsection{Optimality Criteria}
In the next step we choose the objective function also called optimality criterion.
Optimization of the different objectives leads to different experimental design outcomes.
The most popular criteria are~\cite{balsa-cantoe.bangaj.r.COMPUTINGOPTIMALDYNAMIC2008, dette_designing_1997, walter_qualitative_1990}:
\begin{itemize}
\item \textbf{D-optimality criterion} maximizes the determinant $\det (F)$ of the \ac{fim}.
Translated to the parameter space, it means that the volume of the confidence region (see Fig.~\ref{Fig2}) (or the geometric mean of all errors) is minimized.
The size of the confidence region is determined by a confidence level, which is typically chosen $90 \%$ or $95 \%$.
The confidence level can be interpreted as a probability value that the best
parameter fit will belong to this area in case of repeating the experiment and estimations multiple times.
The confidence region is usually presented as a confidence ellipsoid.
D-optimality is the most widely used criterion and is suitable even in case of such parameter transformations as rescaling.
\item \textbf{E-optimality criterion} maximizes the smallest eigenvalue $\lambda_{\min}$ of the \ac{fim}, which is the same as reducing only the largest estimation error.
\item \textbf{A-optimality criterion} maximizes the sum of all eigenvalues $\sum_i \lambda_i$, which can be interpreted as minimizing the algebraic mean of all errors.
\item \textbf{Modified E-optimality criterion} maximizes the ratio between the minimal and maximal eigenvalue $\lambda_{\min} / \lambda_{\max}$ and reduces correlation between two parameters corresponding to these eigenvalues.
\end{itemize}
Each of the criteria has its pros and cons, so the reader should have a closer look at the various properties of these criteria, for instance, in Franceschini's and Macchietto's paper~\cite{franceschiniModelbasedDesignExperiments2008}.
The geometrical interpretation of the criteria using the confidence region is shown in Fig.~\ref{Fig2}.
%
%
% ########################################################################
\subsubsection{Optimization}
After defining the objective function based on the \ac{fim} the next step is to optimize the chosen optimality criteria.
Finding the experimental design corresponds to finding the global maximum of the objective, which can be formulated as an optimal control problem~\cite{espie_optimal_1989}.
This problem has been widely studied, and multiple numerical solution algorithms for both local and global optimization were introduced~\cite{esposito_global_2000, bangaImprovingFoodProcessing2003, ali_numerical_1997, runarsson_stochastic_2000}.
In the supplied toolbox {\it eDPM}, three methods were implemented: differential evolution, basin-hopping, and the so-called “brute force” method.
We give in the following a brief summary of the implemented methods.
For more information the reader should turn to the available literature, e.g.,~\cite{stornDifferentialEvolutionSimple1997, wales_global_1997}.\\
%
The differential evolution algorithm developed by Storn and Price (1996)~\cite{stornDifferentialEvolutionSimple1997} is one of the stochastic global optimization methods appropriate for nonlinear dynamic problems.
Such types of algorithms have mild computational load but one cannot be sure that the absolute optimum is reached~\cite{balsa-cantoe.bangaj.r.COMPUTINGOPTIMALDYNAMIC2008}.
This method is rather simple and straightforward, does not require any gradient calculation and is easy to parallelize.
Another advantage is its fast convergence and robustness at numerical optimization~\cite{babu_differential_2007} (see \textbf{Note~\ref{note-differential-evolution-details}}).\\
%
The second implemented method is basin-hopping developed by David Wales and Jonathan Doye~\cite{wales_global_1997}.
This iterative algorithm combines Monte-Carlo and local optimization (see \textbf{Note~\ref{note-basin-hopping-details}}).
For a deeper understanding of the Monte-Carlo minimization, the reader can refer to the papers~\cite{li_monte_1987, beichl_metropolis_2000}.\\
%
And lastly, a simple brute force method was implemented as well.
It is a grid search algorithm calculating the objective function value at each point of a multidimensional grid in a chosen region.
This method is suitable for discrete optimization with limited number of grid values.
However, the downside of this technique is its slowness, inefficiency and long computational times, noticeable for higher numbers of possible discrete solutions~\cite{scipybrute}.
The papers of Banga \etal~\cite{banga_global_1996, banga_dynamic_2005, bangaImprovingFoodProcessing2003} list additional information about this topic.
%
% ########################################################################
\subsubsection{Optimization Code}
To run the optimization procedure the following function can be used:
\begin{minted}[linenos=false]{python}
fsr = find_optimal(fsm)
\end{minted}
This code will produce the output shown in Fig.~\ref{Fig4} and \ref{Fig5}.
The user can choose between any of the three previously explained optimization procedures.
Any optimization argument of the aforementioned routines can be specified.
The default arguments of the \mintinline[bgcolor=white,style=emacs]{python}{"scipy_differential_evolution"} optimization method are given by the scipy default arguments~\cite{scipydiffev}.
For our specific use-case, we advise to experiment with the options presented in Code Sample~\ref{code:fsr}.
\begin{code}[h]
\begin{minted}[firstnumber=96]{python}
fsr = find_optimal(
# Required argument: The model to optimize
fsm,
# The optimization method
optimization_strategy="scipy_differential_evolution",
# Our custom options
criterion=fisher_determinant,
relative_sensitivities=True,
# Options from scipy.optimize.differential_evolution
recombination=0.7,
mutation=(0.1, 0.8),
workers=-1,
popsize=10,
polish=False,
)
\end{minted}
\caption{Define the optimization conditions and calculate the resulting experimental design.}
\label{code:fsr}
\end{code}
A full list of optional arguments can be seen in the \href{https://docs.scipy.org/doc/scipy/reference/optimize.html#global-optimization}{scipy documentation}~\cite{virtanenSciPyFundamentalAlgorithms2020}.
In addition, there are some interesting optimization options such as the optional arguments \mintinline[bgcolor=white,style=emacs]{python}{relative_sensitivities,criterion}, which are responsible for using relative sensitivities and specifying the optimality criterion as explained in the previous section.
Please view the \href{https://spatial-systems-biology-freiburg.github.io/FisInMa/}{full documentation} for explanation.
The resulting class \mintinline[bgcolor=white,style=emacs]{python}{fsr} contains all definitions, optimal experimental setup (see \textbf{Note~\ref{note-timepoints}}) and information on the optimization process.
For the parameter estimation the initial values $\mbx(t_0)$ of the system are required.
%
% ########################################################################
\subsubsection{Identifiability}
Before using the \acl{oed} generated in the preceding section, the reader needs to check if the parameters of the system are identifiable, i.e., to examine if it is possible to obtain a unique solution for the parameters from the optimization.
It can happen that a subset of the parameters is non-identifiable.
The non-identifiability can be due to the model structure or observables (structural non-identifiability) or insufficient data (practical non-identifiability)~\cite{guillaume_introductory_2019, wieland_structural_2021, walter_identifiability_1996}.
Structural non-identifiability should be avoided as it results in at least one parameter which cannot be determined.
Fortunately, there is a quick and easy way to check whether the system is structural non-identifiable by calculating the rank of the sensitivity matrix~\cite{miao_identifiability_2011, stigter_fast_2015}.
For an identifiable system, the rank should coincide with the number of estimated parameters.
This method allows for the analysis of the local structural identifiability for already chosen inputs and times.
Hence, we can check if the final result of the optimization routine is able to identify every parameter by this method.
This does not replace a systematic inspection of the model but at least we can validate our results.
If a non-identifiability is detected, the structure of the model or the number of measurements does not allow determination of the involved parameters.
This means we need to adjust our description and check all supplied parameters for their validity.
On the other hand, there are also {\it a priori} methods, e.g., using Lie group theory, that only require the model definition~\cite{wieland_structural_2021}.
However, due to the extensiveness of the topic we constrain ourselves to the method discussed above.
Practical non-identifiability results in large confidence intervals~\cite{holmberg_practical_1982, miao_identifiability_2011, kreutzProfileLikelihood2013}.
Methods such as profile-likelihood~\cite{wieland_structural_2021} are able to test for these cases but they require experimental data.
%
% ########################################################################
\subsubsection{Identifiability Code}
Using the {\it eDPM} package, the reader can quite easily check if the resulting experimental design is structurally identifiable.
For this, one can call the function that compares the rank of the sensitivity matrix with the number of parameters as in Code Sample~\ref{code:identifiability_check}.
\begin{code}[h]
\begin{minted}[firstnumber=112]{python}
# Check the structural identifiability
check_if_identifiable(fsr)
\end{minted}
\caption{Check the identifiability of the optimization result.}
\label{code:identifiability_check}
\end{code}
The function returns \mintinline[bgcolor=white,style=emacs]{python}{True}, if the identifiability condition is reached, and \mintinline[bgcolor=white,style=emacs]{python}{False} otherwise.
If this test is not passed, the reader should either reconsider the model structure including the definition of parameters or increase the number of numerically optimized measurement conditions (inputs or times).
From a mathematical point of view, we need at least as many measurement points as we have parameters in our system (see \textbf{Note~\ref{note-structural-identifiability}}).
%
% ########################################################################
\subsubsection{Plotting, Json, etc.}
Our package also provides the option to save results into a Json file and automatically plot results for the ode solutions, observables and sensitivities (see Code Sample~\ref{code:plot_save}).
\begin{code}[h]
\begin{minted}[firstnumber=115]{python}
plot_all_observables(fsr)
plot_all_sensitivities(fsr)
json_dump(fsr, "baranyi.json").
\end{minted}
\caption{Plot the observables and sensitivities calculated for the \ac{oed} and save the result in Json format.}
\label{code:plot_save}
\end{code}
%
%
%
% ########################################################################
% ########################################################################
\subsection{Examples}
\subsubsection{Baranyi-Roberts Model - Single Species}
As a summary of this tutorial, we would like to present the resulting output of our package.
To this end we study the growth of a bacterial colony consisting of a single species and describe it mathematically using the Baranyi-Roberts model given by Eqs.~\ref{eq:ode_BaranyiRoberts},\ref{eq:RatkowskyModel}.
We aim to estimate the parameters of the model and determine optimal measurement conditions.
The observable is the bacterial count, i.e., the first component of the state variable vector $y = x_1$.
As an additional constraint we need to consider that the available climate chambers can only operate the temperature in the range from 2 to 12 degrees, and the experiment should not take longer than 100 hours.
The parameter values and initial values of the system are taken in accordance with the literature~\cite{gospavic_mathematical_2008} and presented in Table~\ref{Table2}.
%
To model realistic uncertainties, we choose the error model given by Eq.~\ref{eq:error_model} with $\gamma_\text{abs}=0.3$ and $\gamma_\text{rel}=0.1$ for all the measurement points.
We identify the three parameters of our system by choosing 4 independent measurement conditions consisting of two temperature-points and two time-points.
In our experimental design optimization routine, we use the D-optimality criterion, relative sensitivities, and the differential evolution optimization method.
The resulting \ac{oed} for the described system is presented in Fig.~\ref{Fig3}, \ref{Fig7}.
The identifiability test confirms the validity of our model setup (Fig.~\ref{Fig6}).
%
%
For the system described above the minimal requirements for passing the identifiability test is using a setup with at least two temperatures and two measurement times per temperature.
The choice of minimal number of measurement times can also be justified when looking at Fig.~\ref{Fig8}.
Here we see that one measurement time is not enough as the determinant is numerically very close to zero indicating a non-identifiability of this system configuration.
%
%
Fig.~\ref{Fig7} shows that each of the chosen optimal time points corresponds to the maximum of at least one local sensitivity curve.
This maximizes the determinant of the \acl{fim}.
%
%
%
\subsubsection{Two-Species Resource Competition}
In this example we discuss a slightly more complicated system and extend the previous example of bacterial growth to a system where two species are present.
The species interact by competitive inhibition because they depend on a common nutrient resource.
Analogously to the previous example, we denote the concentration of the two different species as $x_1$ and $y_1$, respectively.
The system can be described by the following set of \acp{ode}:
\begin{equation}
\begin{cases}
\dot x_1 = \alpha_x R x_1 \\
\dot y_1 = \alpha_y R y_1 \\
\dot R = -\frac{R}{n^\text{max}}(\alpha_x x_1+\alpha_y y_1),
\label{eq:Lotka_Volterra_2species_1}
\end{cases}
\end{equation}
where $\alpha_x, \alpha_y$ are the time and temperature dependent growth rates for population $x_1$, $y_1$, and $R$ represents the nutrient pool concentration, respectively.
Using the conservation quantity $x_1 + y_1 + n^\text{max}R = n^\text{max}$ this system can be reduced to:
\begin{equation}
\begin{cases}
\dot x_1 = \alpha_x x_1\left(1-\frac{x_1+y_1}{n^\text{max}}\right) \\
\dot y_1 = \alpha_y y_1\left(1-\frac{x_1+y_1}{n^\text{max}}\right).
\end{cases}
\label{eq:Lotka_Volterra_2species_2}
\end{equation}
Based on the one species case for the growth rates we use~\cite{baranyiDynamicApproach1994,ratkowsky_relationship_1982}
\begin{eqnarray}
\alpha_x (t, T) = b_x^2 (T - T^\text{min}_x)^2 \frac{x_2(t)}{x_2(t) + 1}\\
\alpha_y (t, T) = b_y^2 (T - T^\text{min}_y)^2 \frac{y_2(t)}{y_2(t) + 1}.
\label{eq:growth_rate_2species}
\end{eqnarray}
Combining Eqs.~\ref{eq:Lotka_Volterra_2species_2} and \ref{eq:growth_rate_2species}, we get a system of four differential equations:
\begin{equation}
\begin{cases}
\dot x_1 = b_x^2 (T - T^\text{min}_x)^2 \frac{x_2}{x_2 + 1} x_1 (1 - \frac{x_1 + y_1}{n^\text{max}})\\
\dot x_2 = b_x^2 (T - T^\text{min}_x)^2 x_2 \\
\dot y_1 = b_y^2 (T - T^\text{min}_y)^2 \frac{y_2}{y_2 + 1} y_1 (1 - \frac{x_1 + y_1}{n^\text{max}})\\
\dot y_2 = b_y^2 (T - T^\text{min}_y)^2 y_2,
\label{eq:model_2species_resource}
\end{cases}
\end{equation}
where $x_2, y_2$ are the concentration of the quantities determining the critical substance (nutrient) needed for growth of the species $x_1$ and $y_1$, respectively.
The values of the parameter vector $\mbp = (n^\text{max}, b_x, T^\text{min}_x, b_y, T^\text{min}_y)$ and the vector of the initial values $\mbx(t_0)=(x_1(t_0), x_2(t_0), y_1(t_0), y_2(t_0))$ are presented in Table~\ref{Table3}.
%
In comparison to the previous example, we now choose the count of the two bacteria species $x_1,y_1$ as observables.
Meanwhile, other optimization arguments were taken the same as in the previous example i.e., we optimize measurement times and temperatures.
The new \ac{oed} is shown in Fig.~\ref{Fig9}.
%
%
The presented design shows that the least needed number of experiments is two with temperatures of $2^\circ$C and $12^\circ$C.
In each experiment two observables that correspond to the concentrations of two bacteria types and two time points were measured (see \textbf{Note~\ref{note-temperature-identifiability}}).
%
%
%
\subsection{Conclusion}
We introduced the concepts of experimental design and identifiability and their practical applications to systems described by \acp{ode}.
In addition, we implemented the numerical calculations in simple python code which we provide within our software {\it eDPM}.
Our toolbox provides a simple way to apply methods of mathematical statistics and optimization.
We demonstrated how to use these methods by applying them to the well-known Baranyi-Roberts model for microbial growth with one and two species.
These results show us how parameters of the model can now be estimated most precisely by designing experiments around the predicted measurement conditions.
In summary, model-based experimental design can accelerate and simplify the planning of efficient experiments for parameter estimation in microbiology.
\section{Notes}
\begin{enumerate}
\item \label{note-order-arguments-python} When using the syntax as shown in Code Sample~\ref{code:model_definition}, the order of arguments does not matter. However, when only using \mintinline[bgcolor=white,style=emacs]{python}{FisherModel(x0, 0.0, ...)}, please pay attention to the order of arguments.
\item \label{note-differential-evolution-details} In differential evolution optimization algorithm an initial population of candidate vectors (sampling times and inputs) is randomly chosen from the region of available values. Then each vector mutates by mixing with one of other candidate vectors: To a chosen vector from the initial population $D_0$, we add a weighted difference between two other randomly chosen vectors from the same set $(D_\text{rand1} - D_\text{rand2})$. A new mutated vector $D_m$ is obtained. The next step is to construct a new trial solution. This is done by choosing the elements of the trial vector either from the initial $D_0$ or the mutated $D_m$ vectors. For each new element of the trial vector a random number is drawn uniformly from the interval $[0, 1)$ and compared to the so-called recombination constant. If this random number is less than the recombination constant, then the trial vector element is chosen from the vector $D_m$, otherwise from $D_0$. The degree of mutation can be controlled by changing the recombination constant; the larger this constant is, the more often vector elements are chosen from $D_m$. Subsequently, the objective function is calculated using the trial vector and the result is compared to result obtained using the initial solution $D_0$; and the best of them is chosen for the next generation. This procedure is repeated for every solution candidate of the initial population, by means of which the new generation is build~\cite{scipydiffev}. The process of population mutation is repeated till a stopping criterion is reached, e.g., the maximum number of generations (steps) is reached or the standard deviation of the candidate vectors is below a certain threshold~\cite{Zielinski_DE}.
\item \label{note-basin-hopping-details} During the iteration step the design vector, i.e., vector of measurement times and inputs, is subject to a random perturbation and then to local minimization. After this, the step is either accepted or rejected. As in a standard Monte-Carlo method, the decision is made using the Metropolis criterion for the objective function~\cite{scipybashop}.
\item \label{note-timepoints} Our toolbox optimizes and provides the optimal experimental timepoints with an assumption that the initial one $t_0$ is always measured. Hence, the reader should keep that in mind and include the initial time $t_0$ as a measurement point in the final \ac{oed}.
\item \label{note-structural-identifiability} The provided identifiability test only allows to exclude structural identifiability but to obtain reasonable confidence intervals for parameter estimates the practical identifiability should be considered as well. Therefore, in reality more experiments are needed to increase the accuracy of the model. Thus, we suggest that the reader considers the optimization result not as a finished design but as a reference pointing to the minimal requirements and the most crucial conditions (inputs and times) for the parameter estimations.
\item \label{note-temperature-identifiability} For both temperatures we require at least two time-points for identifiability. Thus, we can estimate the parameters by performing experiments according to this design.
\end{enumerate}
%
%
%
\section*{Figures}
\begin{figure}[H]
\centering
\includegraphics[scale=0.3]{Figures/Fig1.pdf}
\caption{The workflow of the iterative process for model-based experimental design for parameter estimation.}
\label{Fig1}
\end{figure}
%
%
\begin{figure}[H]
\centering
\includegraphics[scale=0.6]{Figures/Fig2.pdf}
\caption{The confidence ellipsoid projection to the ($p_1$, $p_2$) parameter space.
The ellipsoid shows the geometrical meaning of the optimality criteria.
The center of the ellipsoid represents the estimated parameter values.
The radii are the uncertainties of the estimates associated with different eigenvalues of the \ac{fim} $\lambda_1$, $\lambda_2$ ($\lambda_1 < \lambda_2$).
The D-optimality aims to minimize the volume of the ellipsoid, the E-optimality minimizes the largest radius, A-optimality minimizes the perimeter of the rectangle that encloses the ellipse (dashed gray line), and, finally, modified E-optimality tends to make the ellipse as spherical as possible.}
\label{Fig2}
\end{figure}
%
%
\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{Experiment 1: $T = 3.5^\circ$C, $t=0 \text{h}, 100$h}
\includegraphics[scale=0.255]{Figures/Fig3a.pdf}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{Experiment 2: $T = 12^\circ$C, $t=0\text{h}, 43\text{h}, 100$h}
\includegraphics[scale=0.255]{Figures/Fig3b.pdf}
\end{subfigure}
\caption{The example of the output of the experimental design optimization procedure for the Baranyi-Roberts model~(Eqs.~\ref{eq:ode_BaranyiRoberts} - \ref{eq:observable}).
The line plot presents the model solution for the observable, and the circles determine the time points suggested by experimental design.
The optimal design is proposed for one observable which is the total bacterial count of the species $x_1$ and consists of two time series where samples are stored at temperatures (a) $T_1 = 3.5^\circ$C and (b) $T_2 = 12^\circ$C.
The corresponding measurement times are (a)~$t_{11}=100$h and (b)~$t_{21}=43$h, $t_{22}=100$h.
The initial time $t_0=0$ is included in the experimental design by definition.}
\label{Fig3}
\end{figure}
%
%
% ###### CHANGE THE MINTED ENVIRONMENT ######
%
%
\definecolor{MintedFgColor}{HTML}{fcfcfc}
\definecolor{MintedBgColor}{HTML}{2c2e30}
\AtBeginEnvironment{minted}{\color{MintedFgColor}}
\begin{figure}[H]
\centering
%\includegraphics[scale=1.]{Figures/Fig4.pdf}
\begin{minted}[bgcolor=MintedBgColor,linenos=false,fontsize=\scriptsize]{text}
=========================== SUMMARY OF FISHER MODEL ============================
================================ ODE FUNCTIONS =================================
├─ode_fun baranyi_roberts_ode
├─ode_dfdx ode_dfdx
├─ode_dfdp ode_dfdp
├─obs_fun unknown
├─obs_dgdx unknown
└─obs_dgdp unknown
================================ INITIAL GUESS =================================
├─ode_x0 [array([50., 1.])]
├─ode_t0 [0.]
├─times [[ 0. 100.]
[ 0. 100.]]
├─inputs [array([ 2., 12.])]
├─parameters (20000.0, 0.02, -5.5)
├─ode_args None
├─identical_times False
└─covariance CovarianceDefinition(rel=1.1, abs=1.3)
============================= VARIABLE DEFINITIONS =============================
├─ode_x0 None
├─ode_t0 None
├─times VariableDefinition(lb=0.0, ub=100.0, n=2)
├─inputs [VariableDefinition(lb=2.0, ub=12.0, n=2)]
├─parameters (20000.0, 0.02, -5.5)
├─ode_args None
├─identical_times False
└─covariance CovarianceDefinition(rel=1.1, abs=1.3)
=============================== VARIABLE VALUES ================================
├─ode_x0 [array([50., 1.])]
├─ode_t0 [0.]
├─times [[ 0. 100.]
│ [ 0. 100.]]
├─inputs [array([ 2., 12.])]
├─parameters (20000.0, 0.02, -5.5)
├─ode_args None
├─identical_times False
└─covariance CovarianceDefinition(rel=1.1, abs=1.3)
================================ OTHER OPTIONS =================================
└─identical_times False
\end{minted}
\caption{Command-line output of {\it eDPM} of all parameters, inputs, \ac{ode} functions and initial values needed to run the optimization method.}
\label{Fig4}
\end{figure}
\begin{figure}
%\includegraphics[scale=1.]{Figures/Fig5.pdf}
\begin{minted}[bgcolor=MintedBgColor,linenos=false,fontsize=\scriptsize]{text}
========================== STARTING OPTIMIZATION RUN ===========================
differential_evolution step 1: f(x)= -72.0062
differential_evolution step 2: f(x)= -152.503
differential_evolution step 3: f(x)= -152.503
differential_evolution step 4: f(x)= -152.689
differential_evolution step 5: f(x)= -176.106
differential_evolution step 6: f(x)= -210.415
differential_evolution step 7: f(x)= -224.958
differential_evolution step 8: f(x)= -240.512
differential_evolution step 9: f(x)= -240.512
differential_evolution step 10: f(x)= -240.512
differential_evolution step 11: f(x)= -246.616
differential_evolution step 12: f(x)= -246.616
differential_evolution step 13: f(x)= -247.332
differential_evolution step 14: f(x)= -248.422
differential_evolution step 15: f(x)= -248.744
differential_evolution step 16: f(x)= -249.315
differential_evolution step 17: f(x)= -251.202
differential_evolution step 18: f(x)= -251.202
differential_evolution step 19: f(x)= -251.229
differential_evolution step 20: f(x)= -251.229
differential_evolution step 21: f(x)= -251.469
differential_evolution step 22: f(x)= -251.469
differential_evolution step 23: f(x)= -251.8
differential_evolution step 24: f(x)= -252.206
differential_evolution step 25: f(x)= -252.217
differential_evolution step 26: f(x)= -252.217
differential_evolution step 27: f(x)= -252.31
============================== OPTIMIZED RESULTS ===============================
================================== CRITERION ===================================
├─fisher_determinant 252.31000727743015
├─sensitivity matrix [[0.02936982 6.00979198 3.68012032]
│ [0.02936994 6.00979942 3.68012487]
│ [0.19345575 8.42296024 2.64726938]
│ [0.99551519 0.10814337 0.0339886 ]]
├─inverse covariance matrix [[0.82338185 0. 0. 0. ]
│ [0. 0.82338186 0. 0. ]
│ [0. 0. 0.82594678 0. ]
│ [0. 0. 0. 0.82634818]]
============================== INDIVIDUAL RESULTS ==============================
Result_0
├─ode_x0 [50. 1.]
├─ode_t0 0.0
├─times [99.9942128 99.99433358]
├─inputs [3.4817324001093413]
└─parameters (20000.0, 0.02, -5.5)
Result_1
├─ode_x0 [50. 1.]
├─ode_t0 0.0
├─times [42.95660959 98.6753861 ]
├─inputs [11.999647652116655]
└─parameters (20000.0, 0.02, -5.5)
======================== DISCRETIZATION PENALTY SUMMARY ========================
├─penalty 1.0
├─penalty_ode_t0 1.0
├─penalty_inputs 1.0
├─penalty_times 1.0
└─penalty_summary {'ode_t0': [], 'inputs': [], 'times': []}
\end{minted}
\caption{Command-line output which summarizes the optimization routine and its results.}
\label{Fig5}
\end{figure}
%
\begin{figure}
%\includegraphics[scale=1.]{Figures/Fig6.pdf}
\begin{minted}[bgcolor=MintedBgColor,linenos=false,fontsize=\scriptsize]{text}
=============================== ANALYSIS RESULTS ===============================
└─structural identifiability True
\end{minted}
\caption{Command-line output showing the result of the identifiability check for the experimental design provided in Fig.~\ref{Fig3}.
The output shows that the experimental design gives a structurally identifiable result and allows parameter estimation of the system.}
\label{Fig6}
\end{figure}
%
%
\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{{\footnotesize $T_1=3.5^\circ$C, parameter $x^\text{max}$}}
\includegraphics[scale=0.25]{Figures/Fig7a.pdf}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{{\footnotesize $T_2=12^\circ$C, parameter $x^\text{max}$}}
\includegraphics[scale=0.25]{Figures/Fig7b.pdf}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{{\footnotesize $T_1=3.5^\circ$C, parameter $b$}}
\includegraphics[scale=0.25]{Figures/Fig7c.pdf}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\centering
\subcaption{{\footnotesize $T_2=12^\circ$C, parameter $b$}}