Optimality for models given by a system of ordinary differential
equation

Last update: 2014-01-10

Authors: Juan M. Rodríguez-Díaz and Guillermo Sánchez-León

Many processes are given by a system of ordinary differential equations, very often without an analytical solution. When there are unknown parameters, that need to be estimated, optimum experimental design approach offers quality estimators for the different objectives of the practitioners. But almost every optimality criteria needs to deal with the linearized model for computing optimal designs, and this can be a great problem when it is not possible to obtain the analytical form of the model. In this work, a procedure for findingoptimal designs for models given as solutions to a system of ordinary differential equations is described. Some important models like the compartmental one, are studied through actual case studies, obtaining the corresponding optimal designs.

$Version

Overview of  compartmental and biokinectic equivalents models

Compartmental analysis has applications in clinical medicine, pharmacokinetics, internal dosimetry, nuclear medicine, ecosystem studies and chemical reaction kinetics.  It can be described as the analysis of a system in terms of compartments which separate the system into a finite number of component parts which are called compartments. Compartments interact through the exchange of species. Species may be a chemical substance, hormone, individuals in a population and so on.  A compartmental system is usually represented by a flow diagram or a block diagram. A general introduction to this theory can be found in Anderson (1983), Godfrey (1983) and Jazquez (1985).

We adopt the convention of  representing compartments with circles or rectangles. The flow into or out of the compartments is represented by arrows.  The (i^th) compartment  of a system of n compartments is labelled i and the size (amount or content) of the component in compartment  i as (Subscript[x, i])(t). The exchange between compartments, or between a compartment and the environment is labeled (Subscript[k, ij]), where i represents the flow from i to j. The environment is usually represented by "0" (zero), so (Subscript[k, i0]) is the fractional excretion coefficient from the i-th compartment to the outside environment.  The  input from the environment into the j-th compartment is called (Subscript[b, j])(t). Environment represents the processes that are outside the system. With regards to the environment, we only need to know the flow, (Subscript[b, j])(t),  into the system from the outside. The (Subscript[k, ij]) are called fractional transfer rate coefficients and they may be a function of different variables or constants.

Figure  represents  the general tricompartmental model with one input and output in each compartment. If we suppose that the substance introduced into the system is a radioactive isotope we mush considered the radioactive decay, it is given by a constant rate represented by λ (decay constant), it is specific for each isotope (obviously λ = 0 if no radioactive substances are presents or they are long life isotopes) . The decay constant can be interpreted by an equal flow going out of the system in each compartment. The model can be represented by the follow ODE System

(Subscript[x, 1])'(t) = (Subscript[b, 1])(t) - ((Subscript[k, 10]) + (Subscript[k, 12]) + (Subscript[k, 13])+λ) (Subscript[x, 1])(t)  + (Subscript[k, 21]) (Subscript[x, 2])(t) + (Subscript[k, 31]) (Subscript[x, 3])(t) 
(Subscript[x, 2])'(t) = (Subscript[b, 2])(t) + (Subscript[k, 12]) (Subscript[x, 1])(t)  - ((Subscript[k, 20]) + (Subscript[k, 21]) + (Subscript[k, 23])+λ) (Subscript[x, 2])(t) + (Subscript[k, 32]) (Subscript[x, 3])(t)
(Subscript[x, 3])'(t) = (Subscript[b, 3]) (t) + (Subscript[k, 13]) (Subscript[x, 1])(t)  + (Subscript[k, 23]) (Subscript[x, 2])(t) - ((Subscript[k, 30]) + (Subscript[k, 31]) + (Subscript[k, 32])+λ) (Subscript[x, 3])(t)

in matrix notation

(({{(Subscript[x, 1])'(t)}, {(Subscript[x, 2])'(t)}, {(Subscript[x, 2])'(t)}}))= (({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), (Subscript[a, 33])}}))(({{(Subscript[x, 1])(t)}, {(Subscript[x, 2])(t)}, {(Subscript[x, 2])(t)}}))+(({{(Subscript[b, 1])(t)}, {(Subscript[b, 2])(t)}, {(Subscript[b, 3])(t)}}))

with

(Subscript[a, 11])=- ((Subscript[k, 10]) + (Subscript[k, 12]) + (Subscript[k, 13])+ λ) ;(Subscript[a, 12]) = (Subscript[k, 21]);(Subscript[a, 13])= (Subscript[k, 31]);
(Subscript[a, 21])=(Subscript[k, 12]);(Subscript[a, 22])=- ((Subscript[k, 20]) + (Subscript[k, 21]) + (Subscript[k, 23])+ λ)  ; (Subscript[a, 23])=(Subscript[k, 32]);
(Subscript[a, 31])=(Subscript[k, 13]); (Subscript[a, 32])= (Subscript[k, 23]); (Subscript[a, 33])=- ((Subscript[k, 30]) + (Subscript[k, 31]) + (Subscript[k, 32])+ λ);

In some circumstances, instead of compartments we have parts o process than can be represented by transference between Different parts, for example, when blood flow or other physiological parameter is measured. In those circumstances the transfer rate constants (Subscript[k, ij]) are associated with physiologically meaningful values that correspond to the measured physiological parameter or may be a function of them. However, if the physiological parameters are constants models are mathematically equivalent with compartmental models. We will refer in general to kinetic models.

The patterns that we have seen can be expanded to systems of n compartments or n state variables (in the case of physiological models). The equation for any compartment i in notation matrix is given by

(OverDot[x])(t) =A x(t) + b(t), t≥0
x(0) = (Subscript[x, 0])

where:

x(t)= (({(Subscript[x, 1])(t), (Subscript[x, 2])(t), ...,(Subscript[x, n])(t)})^T) is a column vector and  (Subscript[x, i])(t) denote the amount or content of species in compartment i at time t.

A is a n×n is usually known as the system matrix given by coefficients (Subscript[a, ij])= g((Subscript[k, ij])) where (Subscript[k, ij]) are constant specifics of each model (in compartmental models they are obtained in the form that has been descried)

b(t) = (({(Subscript[b, 1])(t), (Subscript[b, 2])(t),...,(Subscript[b, n])(t)})^T) is a column vector where {(Subscript[b, i])(t)} is the input rate into compartment i from outside system.

x(0)  =  {(Subscript[x, 1])(0), (Subscript[x, 2])(0),  ..., (Subscript[x, n])(0)}Tare the initial conditions, so (Subscript[x, i])(0), represents the amount or content of species in compartment i at time  t = 0.

The solution of eq() when the coeffs (Subscript[a, ij]) are constants is eq()

x(t) =(Subscript[x, 0]) Exp(A t) +Exp(A t) * b(t)

where * denotes convolution,

Exp(A t) * b(t) = (∫_0^t)Exp[(t-τ)A] b(τ)dτ

Sometime Laplace transforms are used to solve eq ().

X(s) =(((s I -A) )^(-1))(Subscript[x, 0])+(((s I -A) )^(-1)) B(s)

where X(s) and B(s) are the Laplace transforms of x(t) and b(t). Then evaluating the inverse Laplace transformation is obtained:

x(t) =(ℒ^(-1))((((s I -A) )^(-1))(Underscript[(Subscript[x, 0]),  ]))+(ℒ^(-1))((((s I -A) )^(-1)) (Underscript[B,  ])(s))

Both methods are used for BIOKMOD software for solving eq.(1) developped for ones of the authors.

In many circumstance some parameters of  biokinetic models represented by eq. (1) (one o more (Subscript[k, ij]) coefficients), that we will call β ={(Subscript[β, 1]),...,(Subscript[β, p]) } , are unknown an they are estimated by fitting experimental data. It means that we measure (Subscript[y, i])(t), being (Subscript[y, i])(t) = (Subscript[x, i])(t, β)+ϵ, for different moments of t: (Subscript[(OverHat[y]), i])((Subscript[t, 0])) ,  ..., (Subscript[(OverHat[y]), i])((Subscript[t, n])). Then we estimated β ={(Subscript[β, 1]),...,(Subscript[β, p]) }  by fitting  (Subscript[y, i])(t) and (Subscript[(OverHat[y]), i])(t).

We propose chose the best moments {(Subscript[t, 0]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} to take the experimental data. It  it can be done using techniques of Optimal Design of Experiment (ODE), in particular we chose apply a  D-optimal design.

In a D-optimal design  the values of {(Subscript[t, 1]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} ((Subscript[t, i]) is the time when the i-th sample should be taken) are given in the points that leads the determinant of the Fisher information matrix (M) to a maximum.  The process to obtain M will be described later (it can be found in Hill, P. D. H., 1980, D-optimal designs for partially nonlinear regression models. Technometrics 22:275–276).

Derivatives

The  standard method to obtain the maximum of Det|M| requires to know the analytical expression of (Subscript[x, i])(t, β). In this paper we propose a method to compute the D-optimal design in biokinetic system descried by ODE with the pattern of eq (2) when analytical expression of (Subscript[x, i])(t, β) cannot be obtained.

In the D-optimal design method the derivatives with respect to β must be evaluated

We used a subscript in parentheses to denote differentiation with respect to a parameter, so ∂x/∂β = (Subscript[x, (β)]).

When in eq(1) A depends of β, but not b(t) and (Subscript[x, 0]) we can get the derivatives of  (Subscript[x, i])(t, β) by differentiation of eq.(1)

(Subscript[(OverDot[x]), (p)])(t) =A (Subscript[x, (p)])(t) +(Subscript[A, (p)]) x(t)

Then, using  eq()

(Subscript[x, (p)])(t) =Exp(A t)  (Subscript[x, (p)])(0) +Exp(A t)*[(Subscript[A, (p)])  x(t) ] =Exp(A t)  (Subscript[x, (p)])(0) +Exp(A t) *[(Subscript[A, (p)])  Exp(A t) (Subscript[x, 0])]+Exp(A t) *(Subscript[A, (p)])Exp(A t) * b(t)

(Subscript[x, (p)])(t) =Exp(A t)  (Subscript[x, (p)])(0) +(Subscript[A, (p)])(Subscript[x, 0])Exp(A t) * Exp(A t)+(Subscript[A, (p)])Exp(A t) *Exp(A t) * b(t)

Comentario: (Subscript[x, (a)])(0) = d(Subscript[x, 0])/da. Así se indica en  Bates, en el parrafo que sigue a la ec (5.14). 
Para calcular X(p)(0) lo que hace es d(Xo)/dp, es decir deriva X(0)= Xo respecto de a. Lo anterior lo corrobora en el ejemplo de la pag. 179 Tetracycline 5.La conclusión es que siempre se debe tomar  X(p)(0)=0, excepto cuando las condiciones iniciales son funcion de p, esto es X(0)= Xo(p).
Para mi tiene sentido fisiológico, si p es un factor de trasferencia no experimentara ningun efecto sobre las condiciones iniciales y por tanto suderivada sera cero. Sin embargo si p es un parametro que afecta a las condiciones iniciales su derivada sera distinto de cero. Por ejemplo: en el modelo pulmonar la deposicion en el pulmon en t = 0 de pendende del amad p. En este caso x(0) si depende del AMAD p y por tanto su valor es distinto de cero y su derivada en t=0 tambien lo es. De hecho en mi opinión no tiene sentido tener que calcular x(p) (t), entonces derivar X(p)(t)
respecto de t y despues calcular el valor que toma en t=0 pues para ello necesitas conocer previamente X(p)(t) que es justo lo que se pretende calcular

As x(0) = (Subscript[x, 0 ]) ⇒ (Subscript[x, (p)])(0)=0 therefore:

(Subscript[x, (p)])(t) =(Subscript[A, (p)])(Subscript[x, 0])Exp(A t) * Exp(A t)+(Subscript[A, (p)])Exp(A t) *Exp(A t) * b(t)

In the the particular case that the input happen in t=0, that is  x(0) = (Subscript[x, 0]) with  b(t)=0 for t>0. Then

(Subscript[x, (p)])(t) =(Subscript[A, (p)])(Subscript[x, 0])Exp(A t) * Exp(A t)

Exp(A t) * Exp(A t)  = (∫_0^t)Exp[(t-τ)A] Exp(A τ) dτ

The iodine model

The model

Let's consider the iodine biokinetic model represented in the figure  (ICRP 78) where compartment 1 is the blood, compartment 2 is the thyroid, compartment 3 is the rest of the body, compartment 4 is the bladder,  3->0, i.e. a transfer from compartment 3 to the environment, represents the output to the gastro intestinal tract (GIT) and 4->0 represents the output, via urine excretion, to the environment. In the follow we won`t consider compartment 4 that is not relevant in our case. We will assume a flow from compartment 1 to outside given by a transfer coefficients (Subscript[k, 10]).  The coefficient transfer values, in (days^(-1)), taken from ICRP 78 are (Subscript[k, 10])= 1.9404,  (Subscript[k, 30])= 0.01155 and (Subscript[k, 31])= 0.0462. We will suppose that (Subscript[k, 12]) and (Subscript[k, 23]) are unknown, although we know that their values will be about  (Subscript[k, 12])=0.8 and (Subscript[k, 23])=0.0078. We wish estimate them taken experiment  data from compartment 1 . The problem consist on decide by DOE the best moment to taken the sample.

Here it is downloaded the package Sysmodel (included in the Biokmod Tool available in http://www3.enusa.es/webMathematica/Public/Docs/biokmod.zip)

In this example it is assumed an input into (Subscript[b, 1])(t) compartment 1 is given by

B={-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.020 (E^(-0.147 t))+0.0194 (E^(-0.093 t)),0,0};

Plot[B[[1]], {t, 0, 3}, AxesLabel->{"days", "r"}]

This kind of input happens in real situations when there is an input from the GIT (Gastro Intestinal) to the blood, for instance if the iodine is intaken by orally. Then (Subscript[b, 1])(t) represents the input from GIT to blood. (When the input happens by inhalation the flow to the blood can be represented by 34.4` (e^(-200.` t))+1.09` (e^(-110.` t))+0.808` (e^(-102.` t))+6.414` (e^(-100.` t))+5.458` (e^(-24.` t))  that is the function used in the article). The initial condition are { 0, 0, 0}.

In this kind of experiment usually is used a isotope of iodine with a decay constant λ (his value depend on the isotope). The coefficients  matrix is:

A= CompartMatrix[3, {{1, 2, k12}, {1, 0, 1.9404}, {2, 3, k23} , {3, 0, 0.01155} , {3, 1, 0.0462}},λ]//Chop//TraditionalForm

ShowODE[{{-1.9404-k12-λ,0,0.0462},{k12,-k23-λ,0},{0,k23,-0.05775-λ}}, {0,0,0}, B, t, x] //TraditionalForm

We represent the evolution of the iodine content in the compartment 1 where the samples will be taken  (using (Subscript[k, 12])=0.8 and (Subscript[k, 23])=0.0078). We will refer to iodine 131 which has a radioactive half-life of 8 days, this meaning that radioactive decay constant λ = ln 2/8.02  (day^(-1)). Then the compartmental matrix is:

iodine131matrix= CompartMatrix[3, {{1, 2, 0.8}, {1, 0, 1.9404}, {2, 3, 0.0078} , {3, 0, 0.01155} , {3, 1, 0.0462}}, Log[2]/8.02]

{ x1[t_],x2[t_],x3[t_]} = {(Subscript[x, 1])[t],(Subscript[x, 2])[t],(Subscript[x, 3])[t]}/.SystemDSolve[iodine131matrix, {0,0,0}, B, t,t, x]//Chop;

LogPlot[x1[t], {t, 0, 50},PlotRange-> All, AxesLabel->{"days", "r"}]

The form of the first interval it is explained because in during this interval   (about 3 days) is happening a input b1(t) to compartment 1 and negligible flow from other compartments reach the compartment 1. In the second interval the opposite happens.

We wish obtain (Subscript[(Subscript[x, i]), (β)]) =  {∂(Subscript[x, i])/∂(Subscript[k, 12]), ∂(Subscript[x, i])/∂(Subscript[k, 23])} , i = {1,2,3} to be used later for computing the Optimal Design.

We are going to apply different methods to obtain (Subscript[x, (β)]) in order to make a comparison.  Each method is computated in a new Mathematica session in order to compare the computation time.

Quit[]

Method 1

Here will be obtain (Subscript[x, (β)]) using eq(5) that we call the method of the extended matrix.

We call:Xa={x1a=((∂(Subscript[x, 1])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),x2a=((∂(Subscript[x, 2])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),x3a=((∂(Subscript[x, 3])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12])))} and Ak12 =(Subscript[A, ((Subscript[k, 12]))]) ;Xb={x1b=((∂(Subscript[x, 1])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23]))),x2a=((∂(Subscript[x, 2])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23]))),x3a=((∂(Subscript[x, 3])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23])))} and Ak23 =(Subscript[A, ((Subscript[k, 23]))])

A= CompartMatrix[3, {{1, 2, k12}, {1, 0, 1.9404}, {2, 3,k23} , {3, 0, 0.01155} , {3, 1, 0.0462}}, Log[2]/8.02] //Chop

B={-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.020 (E^(-0.147 t))+0.0194 (E^(-0.093 t)),0,0};

X= {x1[t], x2[t],x3[t]};

eq1=Thread[D[X,t]==A.X + B]

ic1={x1[0] == 0, x2[0] == 0, x3[0] == 0 };

Xa= {x1a[t], x2a[t],x3a[t]};

Xb= {x1b[t], x2b[t],x3b[t]};

Ak12=D[A, k12]

Ak23 = D[A, k23]

eq2=Thread[D[Xa,t]==Ak12.X + A.Xa]

eq3=Thread[D[Xb,t]==Ak23.X + A.Xb];

ic2={x1a[0] == 0, x2a[0] == 0, x3a[0] == 0 };

ic3={x1b[0] == 0, x2b[0] == 0, x3b[0] == 0 };

Now eq1 and eq2, and the respective initial conditions ic1 and ic2 are combined obtaining the eq4

eq4=Join[eq1, eq2 , ic1, ic2]//Chop

eq5=Join[eq1, eq3 , ic1, ic3]//Chop

Them eq4 and eq5 can be solved  when specific values of t, (Subscript[k, 12]) and (Subscript[k, 23]) are given

fa[a_?NumberQ, b_?NumberQ, t1_?NumberQ]:=Evaluate[{x1a[t],x2a[t],x3a[t]}/.NDSolve[Evaluate[eq4/.{k12->a, k23->b}],Join[X, Xa], {t, 0, 100}]]/.t->t1

fb[a_?NumberQ, b_?NumberQ, t1_?NumberQ]:=Evaluate[{x1b[t],x2b[t],x3b[t]}/.NDSolve[Evaluate[eq5/.{k12->a, k23->b}],Join[X, Xb], {t, 0, 100}]]/.t->t1

Then  (Subscript[x, (p)])(t) :{x1a[t]  x1b[t]}.

X1[a_,b_, ti_]:={ fa[a,b, ti][[1,1]], fb[a,b, ti][[1,1]]}

So for a=k12 = 0.80; b= k23= 0.0078, then (Subscript[x, (p)])(t) :{x1a[t], x2a[t],x3a[t]} for t = {1,10,30} are (The computation time, in s, is the first value of the Output.)

Map[X1[0.80, 0.0078, #]&,{1,10,30}]//AbsoluteTiming

Other option faster than the previous is using the new function (Mathematica 9 or later is requiered)   ParametricNDSolve.

eq4a= ParametricNDSolve[eq4,{x1,x2,x3,x1a,x2a,x3a}, {t, 0, 100},{k12,k23}];

eq5a= ParametricNDSolve[eq5,{x1,x2,x3,x1b,x2b,x3b}, {t, 0, 100},{k12,k23}];

fa[a_?NumberQ, b_?NumberQ, t_?NumberQ]:=x1a[a, b][t]/.eq4a

fb[a_?NumberQ, b_?NumberQ, t_?NumberQ]:=x1b[a, b][t]/.eq5a

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

So for a=k12 = 0.80; b= k23= 0.0078, then (Subscript[x, (p)])(t) :{x1a[t], x2a[t],x3a[t]} for t = {1,10,30} are (The computation time, in s, is the first value of the Output.)

Map[X1[0.80, 0.0078, #]&,{1,10,30}]//AbsoluteTiming

Quit[]

Method 2

This method give the solution as function of the unknown parameters (k12 and k23) then when are given specific values of k12 and k23 the ODE is  numerically solved and the derivatives (Subscript[x, (β)]) evaluated numerically in each point. To solve the ODE is used  the Biomod function SystemDSolve (it applies the Mathematica function NDSolve).

A[k12_, k23_]= CompartMatrix[3, {{1, 2, k12}, {1, 0, 1.9404}, {2, 3,k23} , {3, 0, 0.01155} , {3, 1, 0.0462}}, Log[2]/8.02] //Chop

model[t1_?NumberQ,k12_?NumberQ, k23_?NumberQ] := model[t1,k12, k23] ={(Subscript[x, 1])[t1],(Subscript[x, 2])[t1],(Subscript[x, 3])[t1]} /.SystemNDSolve[A[k12,k23],{0, 0, 0}, {-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.020 (E^(-0.147 t))+0.0194 (E^(-0.093 t)),0,0} ,{t, 0, 100},t1,x];

Note: Instead of the before function can be used the following function that solved analytically the ODE when k12 and k23 take numeric values but the take of computation is too longer

model[t1_?NumberQ,k12_?NumberQ, k23_?NumberQ] := model[t1,k12, k23] ={(Subscript[x, 1])[t1],(Subscript[x, 2])[t1],(Subscript[x, 3])[t1]} /.SystemDSolve[{...} ,t,t1,x];

We use the package NumericalCalculus to compute the numerical derivations

Needs["NumericalCalculus`"]

We call fa(t1, (Subscript[k, 12]), (Subscript[k, 23]))={((∂(Subscript[x, 1])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),((∂(Subscript[x, 2])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),((∂(Subscript[x, 3])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12])))} and fb(t1, (Subscript[k, 12]), (Subscript[k, 23]))={((∂(Subscript[x, 1])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23]))),((∂(Subscript[x, 2])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23]))),((∂(Subscript[x, 3])(t,(Subscript[k, 23])))/(∂(Subscript[k, 23])))} (for convenience we write a1, b2 instead of  (Subscript[k, 12]), (Subscript[k, 23]))

Note that sometime the option ND[model[ti,x,b1],x,a1,Scale->.01] should be  used (see ND Help)

fa[a1_,b1_, ti_]:=ND[model[ti,x,b1],x,a1]

fb[a1_,b1_,ti_]:=ND[model[ti,a1,y],y,b1]

Now we test the method using the same values that in the previous example, that is  k12 = 0.80; k23 = 0.0078;

X1[a_,b_,ti_]:={ fa[a,b,ti][[1]], fb[a,b,ti][[1]]}

Map[X1[0.80, 0.0078, #]&,{1,10,30}]//AbsoluteTiming

The solution  is almost the same that the obtained using Method 1 and the computation time is a bit bigger.

Quit[]

Method 3

This method is similar to Method 2 but here is used the new Mathematica (9 o later) funcion ParametricNDSolve.

A= CompartMatrix[3, {{1, 2, k12}, {1, 0, 1.9404}, {2, 3,k23} , {3, 0, 0.01155} , {3, 1, 0.0462}}, Log[2]/8.02] //Chop

B={-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.020 (E^(-0.147 t))+0.0194 (E^(-0.093 t)),0,0};

eqs =ShowODE[A,  {0,0,0}, B, t, x]

sol= ParametricNDSolve[eqs,{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t, 0, 100},{k12,k23}]

fa[a1_?NumberQ, b_?NumberQ, t_?NumberQ]:=D[(Subscript[x, 1])[a,b],a][t]/.a->a1/.sol

fb[a_?NumberQ, b1_?NumberQ, t_?NumberQ]:=D[(Subscript[x, 1])[a,b],b][t]/.b->b1/.sol

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

So for a=k12 = 0.80; b= k23= 0.0078, then (Subscript[x, (p)])(t) :{x1a[t], x2a[t],x3a[t]} for t = {1,10,30} are (The computation time, in s, is the first value of the Output.)

Map[X1[0.80, 0.0078, #]&,{1,10,30}]//AbsoluteTiming

Quit[]

The solution  is almost the same that the obtained using Method 1 and 2 and the computation time is similar.

Method 4

In this case we will obtain again {((∂(Subscript[x, 1])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),((∂(Subscript[x, 2])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12]))),((∂(Subscript[x, 3])(t,(Subscript[k, 12])))/(∂(Subscript[k, 12])))} but in this case we will use  eq(8) with   (Subscript[x, 0])=0

(Subscript[x, (p)])(t) =Exp(A t) *(Subscript[A, (p)])Exp(A t) * b(t)

A={{-2.026827329246876-k12,0,0.0462},{k12,-0.08642732924687598-k23,0},{0,k23,-0.14417732924687598}};

B={-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.020 (E^(-0.147 t))+0.0194 (E^(-0.093 t)),0,0};

func[t1_?NumberQ,a_?NumberQ, b_?NumberQ, k_]:= Module[{m1, AExp},
m1= MatrixExp[A t]/.{k12->a, k23->b}//ExpandAll//Chop;
AExp=Map[Integrate[#1, {tau, 0, t}] &, Evaluate[m1 /. t -> t - tau].Evaluate[D[A,k ].Evaluate[m1 /. t -> tau]], {1}];
Map[Integrate[#1, {tau, 0, t}] &, Evaluate[AExp /. t -> t - tau].Evaluate[B/. t -> tau], {1}]/.t->t1]

The solution of fa is wrong (see the solution obtained with Method 1 and 2) and the time of computation is too long.  We discart this method

{func[30, 0.80, 0.0078, k12],func[30, 0.80, 0.0078, k23]}//AbsoluteTiming

Note that this function is equivalente to the previous function: X1[a_,b_,ti_]:={fa[a,b,ti][[1]], fb[a,b,ti][[1]]}

Conclusion:  Method 1 and 3 are very fast and they are  also  the easiest for programming. We will compare both methods in a OED

Quit[]

Optimal experiment design

We will suppose that (Subscript[k, 12]) and (Subscript[k, 23]) are unknown, although we know that their values will be about  (Subscript[k, 12])=0.8 and (Subscript[k, 23])=0.0078. We wish estimate them taken experiment  data from compartment 1 . The problem consist on decide by DOE the best moment to taken the sample. We will use   D-optimal design.

Method 1

Here we will the optimal design experiment computing the derivatives using the method 1 that we have yet described

We wish  find  t :{(Subscript[t, 0]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} of the model given by eq. (4). (or (5))  using D-optimal design when the analitycal expression of (Subscript[x, 1])(t,  a,b) can not be found. [f(t, β) =(Subscript[x, 1])(t,  (Subscript[k, 12]),(Subscript[k, 23]))]

1.- It is defined a model  f(t, β) where the unknown parameters are β = {a,b}. In our case we call  β = {(Subscript[k, 12]),(Subscript[k, 23])}. [We write eq4a and eq5a obtained when the method 1 has beed described ]

eq4a= ParametricNDSolve[{x1'[t]==-27.13 (E^(-24.08 t))+27.13 (E^(-2.86 t))-0.02 (E^(-0.147 t))+0.0194 (E^(-0.093 t))+(-2.026827-k12) x1[t]+0.0462 x3[t],x2'[t]==k12 x1[t]+(-0.086427-k23) x2[t],x3'[t]==k23 x2[t]-0.1441773 x3[t],x1a'[t]==-x1[t]+(-2.0268-k12) x1a[t]+0.0462` x3a[t],x2a'[t]==x1[t]+k12 x1a[t]+(-0.086427-k23) x2a[t],x3a'[t]==k23 x2a[t]-0.144177` x3a[t],x1[0]==0,x2[0]==0,x3[0]==0,x1a[0]==0,x2a[0]==0,x3a[0]==0},{x1,x2,x3,x1a,x2a,x3a}, {t, 0, 100},{k12,k23}];

eq5a= ParametricNDSolve[ {x1'[t]==-27.13 (E^(-24.08` t))+27.13 (E^(-2.86 t))-0.02 (E^(-0.147 t))+0.0194 (E^(-0.093 t))+(-2.0268-k12) x1[t]+0.0462 x3[t],x2'[t]==k12 x1[t]+(-0.0864-k23) x2[t],x3'[t]==k23 x2[t]-0.1442 x3[t],x1b'[t]==(-2.0268-k12) x1b[t]+0.0462` x3b[t],x2b'[t]==k12 x1b[t]-x2[t]+(-0.0864-k23) x2b[t],x3b'[t]==x2[t]+k23 x2b[t]-0.14418 x3b[t],x1[0]==0,x2[0]==0,x3[0]==0,x1b[0]==0,x2b[0]==0,x3b[0]==0},{x1,x2,x3,x1b,x2b,x3b}, {t, 0, 100},{k12,k23}];

2.-Now it is computed  ∇(f(t), {a, b}) = {((df(t))/da),((df(t))/db)},

fa[a_?NumberQ, b_?NumberQ, t_?NumberQ]:=x1a[a, b][t]/.eq4a

fb[a_?NumberQ, b_?NumberQ, t_?NumberQ]:=x1b[a, b][t]/.eq5a

3.- We need to  define the number of points n to be used in the optimal design.

4.- It is evaluated ∇(f(t), β) at points t:{(Subscript[t, 0]),..., (Subscript[t, n])}, obtaining X = {(Subscript[X, 1]),..., (Subscript[X, p])} with (Subscript[X, 1]) = {((df((Subscript[t, 0])))/(Subscript[dβ, 1])), ...,((df((Subscript[t, n])))/(Subscript[dβ, 1])) }, ..., (Subscript[X, p]) = {((df((Subscript[t, 0])))/(Subscript[dβ, p])), ...,((df((Subscript[t, n])))/(Subscript[dβ, p])) }  . Because the sample will be taken in compartment 1 , we extract of fa and fb the derivatives corresponding to (Subscript[x, 1])(t)

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

5.- A typical election for compute the covariance matrix is assumed that  that the relationship between samples decays exponentially with increasing time-distance between them, that is  Γ = {(Subscript[l, ij])} with (Subscript[l, ij])= exp {ρ|(Subscript[t, j]) -(Subscript[t, j])|}.For computational purpose we have found more appropriate to use the distance (Subscript[d, i]) = (Subscript[t, i])-(Subscript[t, i-1]), instead of  (Subscript[t, i]), then (Subscript[t, i])= (∑_i)(Subscript[d, i]) being (Subscript[d, 0]) =(Subscript[t, 0]) . That is for a two points design . We suppose a 3-points design. The first is defined by the user 

 Γ where

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2)))},{(E^(-ρ d1)),1,(E^(-ρ d2))},{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1}};

6.-  Now it is computed the covariance matrix Σ = (σ^2)  Γ

Σ = (σ^2)*Γ;

We assume

ρ=1;σ =1;

We will also  need give the initial values of  β  the standard deviation of the measures. We also assumed k12= 0.80, k23=0.0078

7.- Then we can obtain the information matrix

M =  (X^T) (Σ^(-1)) X

m := X . Inverse[Σ].  Transpose[X];

m1[ti_] :=Transpose[Map[X1[0.80, 0.0078, #]&,ti] ]. Inverse[Σ].Map[X1[0.80, 0.0078, #]&,ti]

8.- Finally  the determinant of the information matrix is maximized as function of d0, d1 and d2. We constrain the d values to a maximun of t=50 becouse to longer time the concentration will be very low (lower than the detection limit)

obj[d0_?NumericQ,d1_?NumericQ,d2_?NumericQ]:=Det[m1[{d0,d1+d0,d0+d1+d2}]]

sol1=NMaximize[ {obj[d0,d1,d2],0<d0<50,0.02<d1<50, 0.02<d2<50},{d0,d1,d2}]//Timing

Quit[]

Method 3

Here we will the optimal design experiment computing the derivatives using the method 3 that we have yet described.

1.- It is defined a model  f(t, β) where the unknown parameters are β = {a,b}. In our case we call  β = {(Subscript[k, 12]),(Subscript[k, 23])}. [We write the ODE of the system obtained when we described the method 3

sol= ParametricNDSolve[{(Subscript[x, 1])'[t]==-27.13` (E^(-24.08` t))+27.13` (E^(-2.86` t))-0.02` (E^(-0.147` t))+0.0194` (E^(-0.093` t))+(-2.026827`-k12) (Subscript[x, 1])[t]+0.0462` (Subscript[x, 3])[t],(Subscript[x, 2])'[t]==k12 (Subscript[x, 1])[t]+(-0.08643`-k23) (Subscript[x, 2])[t],(Subscript[x, 3])'[t]==k23 (Subscript[x, 2])[t]-0.1442` (Subscript[x, 3])[t],(Subscript[x, 1])[0]==0,(Subscript[x, 2])[0]==0,(Subscript[x, 3])[0]==0},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t, 0, 100},{k12,k23}];

2.-Now it is computed  ∇(f(t), {a, b}) = {((df(t))/da),((df(t))/db)},

fa[a1_?NumberQ, b_?NumberQ, t_?NumberQ]:=D[(Subscript[x, 1])[a,b],a][t]/.a->a1/.sol

fb[a_?NumberQ, b1_?NumberQ, t_?NumberQ]:=D[(Subscript[x, 1])[a,b],b][t]/.b->b1/.sol

3.- We need to  define the number of points n to be used in the optimal design.

4.- It is evaluated ∇(f(t), β) at points t:{(Subscript[t, 0]),..., (Subscript[t, n])}, obtaining X = {(Subscript[X, 1]),..., (Subscript[X, p])} with (Subscript[X, 1]) = {((df((Subscript[t, 0])))/(Subscript[dβ, 1])), ...,((df((Subscript[t, n])))/(Subscript[dβ, 1])) }, ..., (Subscript[X, p]) = {((df((Subscript[t, 0])))/(Subscript[dβ, p])), ...,((df((Subscript[t, n])))/(Subscript[dβ, p])) }  . Because the sample will be taken in compartment 1 , we extract of fa and fb the derivatives corresponding to (Subscript[x, 1])(t)

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

5.- A typical election for compute the covariance matrix is assumed that  that the relationship between samples decays exponentially with increasing time-distance between them, that is  Γ = {(Subscript[l, ij])} with (Subscript[l, ij])= exp {ρ|(Subscript[t, j]) -(Subscript[t, j])|}.For computational purpose we have found more appropriate to use the distance (Subscript[d, i]) = (Subscript[t, i])-(Subscript[t, i-1]), instead of  (Subscript[t, i]), then (Subscript[t, i])= (∑_i)(Subscript[d, i]) being (Subscript[d, 0]) =(Subscript[t, 0]) . That is for a two points design . We suppose a 3-points design. The first is defined by the user 

 Γ where

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2)))},{(E^(-ρ d1)),1,(E^(-ρ d2))},{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1}};

6.-  Now it is computed the covariance matrix Σ = (σ^2)  Γ

Σ = (σ^2)*Γ;

We assume

ρ=1;σ =1;

We will also  need give the initial values of  β  the standard deviation of the measures. We also assumed k12= 0.80, k23=0.0078

7.- Then we can obtain the information matrix

M =  (X^T) (Σ^(-1)) X

m := X . Inverse[Σ].  Transpose[X];

m1[ti_] :=Transpose[Map[X1[0.80, 0.0078, #]&,ti] ]. Inverse[Σ].Map[X1[0.80, 0.0078, #]&,ti]

8.- Finally  the determinant of the information matrix is maximized as function of d0, d1 and d2. We constrain the d values to a maximun of t=50 becouse to longer time the concentration will be very low (lower than the detection limit)

obj[d0_?NumericQ,d1_?NumericQ,d2_?NumericQ]:=Det[m1[{d0,d1+d0,d0+d1+d2}]]

sol1=NMaximize[ {obj[d0,d1,d2],0<d0<50,0.02<d1<50, 0.02<d2<50},{d0,d1,d2}]//Timing

Conclusion:  The time of computation of Method 1 and 3 are practically the same, method 1 a bit faster tham method 3,, but Method 3 is the easiast for programming

Here it is shown graphically the iteration process of d1 and de

d0 = 0.748664;

FindMaximum[{obj[d0,d1,d2],0.02<d1<10, 0.02<d2<10},{{d1,5},{d2,6}},StepMonitor:>Print[{"d1:",d1, "d2:",d2}]]

pts = Reap[FindMaximum[ {obj[d0,d1,d2],0.02<d1<10, 0.02<d2<10},{{d1,5},{d2,6}},StepMonitor:>Sow[{d1,d2}]]][[2,1]];

ContourPlot[obj[d0,d1,d2]//Minus,{d1,0.02,10},{d2, 0.02,10},Epilog->{Red,Line[pts],Point[pts]}]

Quit[]

Biokinetic model of Ciprofloxacin and Ofloxacin

The model

The following differential equations result from the model  ofA. Sánchez-Navarro, C. Casquero, and M. Weiss, 'Distribution of Ciprofloxacin and Ofloxacin in the Isolated Hindlimb of the Rat', Pharmaceutical Research, 16: 587-591 (1999):

(Subscript[c, out])'[t]+((Q/(Subscript[V, p]))+ (PS/(Subscript[V, p])))(Subscript[c, out])[t] -(PS/(Subscript[V, p])) (Subscript[c, Tu])[t]=(Q/(Subscript[V, p])) (Subscript[c, in])[t]
(Subscript[c, Tu])'[t] +((PS/(Subscript[V, Tu])) +(Subscript[k, on]) )(Subscript[c, Tu])[t]-(PS/(Subscript[V, Tu]))(Subscript[c, out])[t]-(Subscript[k, off])((Subscript[V, Tb])/(Subscript[V, Tu]))(Subscript[C, Tb])[t]= 0
(Subscript[c, Tb])'[t]+(Subscript[k, off])(Subscript[c, Tb])[t]- (Subscript[k, on])((Subscript[V, Tu])/(Subscript[V, Tb]))(Subscript[c, Tu])[t]=0

Initial conditions: (Subscript[c, out])[0]= 0,(Subscript[c, Tu])[0]=0, (Subscript[c, Tb])[0]=0

(Subscript[c, i])(t) , with i = {out, Tu, Tb}, represents the concentration of a sustance in differentes part of the models. We consider the case where with Q =  3 mL (min^(-1)), (Subscript[V, Tu])=6.411, (Subscript[V, p]), 0.973, (Subscript[V, Tb])=1,  and P S=2.714, then replacing (Subscript[V, Tu]), (Subscript[V, p]),  (Subscript[V, Tb]), P S, Q for their values where for convenients we call  (Subscript[x, 1])(t), (Subscript[x, 2])(t) and (Subscript[x, 3])(t) intead of (Subscript[c, out])(t), (Subscript[c, Tu])(t)and (Subscript[c, Tb])(t)

({{(Subscript[x, 1])'(t)==-5.872559095580678` (Subscript[x, 1])(t)+2.7893 (Subscript[x, 2])(t) +3.083247687564234` (Subscript[c, int])(t)}, {(Subscript[x, 2])'(t)==0.42333489315239436` (Subscript[x, 1])(t)+(-(Subscript[k, on])-0.42333489315239436`) (Subscript[x, 2])(t)+0.15598 (Subscript[k, off]) (Subscript[x, 3])(t)}, {(Subscript[x, 3])'(t)==6.411` (Subscript[k, on]) (Subscript[x, 2])(t)-(Subscript[k, off]) (Subscript[x, 3])(t)}, {(Subscript[x, 1])(0) =(Subscript[x, 2])(0)==(Subscript[x, 3])(0)==0}})

We use  (Subscript[c, in])(t) = 13610.10776942386` t  (e^(-11.215992508521028` t))  (According to G. Sanchez Biokmod: A Mathematica toolbox for modeling Biokinetic Systems”. Mathematica in Education and Research: 10 (2) 2005. ISSN/ISBN: 1096-3324  using the experimental data obtained by of A. Sánchez-Navarro et al. )

Then, on notation matrix:

x(t) =(Subscript[x, 0]) Exp(A t) +Exp(A t) * b(t)

where

x(t)= (({(Subscript[x, 1])(t), (Subscript[x, 2])(t),(Subscript[x, 3])(t)})^T)

Optimal experiment design

This model (eq9) is represented in Mathematica  by the ODE

eqox={ x1'[t] == -5.87256  x1[t] + 2.78931 x2[t] + 41963.3 t Exp[-11.216 t], x2'[t] == 0.423335 x1[t] + (-kon - 0.423335) x2[t] + 0.155982  koff x3[t],
 x3'[t] == 6.411 kon x2[t] - koff x3[t],
 x1[0] == 0, x2[0] == 0, x3[0] == 0};

We want to estimate the values of (Subscript[k, on]) and (Subscript[k, off]) making an experiment that consist on in measure the concentration of a compount at (Subscript[(OverHat[x]), 1])(t)  for t :{(Subscript[t, 0]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} an then estimate (Subscript[k, on]) and (Subscript[k, off]) by fitting  (Subscript[x, 1])(t,  (Subscript[k, on]),(Subscript[k, off])) .

The problem is can not be found an analitycal expresion of (Subscript[x, 1])(t,  (Subscript[k, on]),(Subscript[k, off])) because the ODE system eq. (4). (or (5)) can not be solved if (Subscript[k, on]),(Subscript[k, off]) are parameters. However  (Subscript[x, 1])(t,  (Subscript[k, on]),(Subscript[k, off])) has solution when  (Subscript[k, on]),(Subscript[k, off]) take numeric values, this fact is used by some nonlinear regresion method  to estimated (Subscript[k, on]),(Subscript[k, off]).

We wish  find  t :{(Subscript[t, 0]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} of the model given by eq. (4). (or (5))  using D-optimal design when the analitycal expresion of (Subscript[x, 1])(t,  (Subscript[k, on]),(Subscript[k, off])) can not be found. [f(t, β) =(Subscript[x, 1])(t,  (Subscript[k, on]),(Subscript[k, off]))]

1.- It is defined a model  f(t, β) where the unkown parameters are β = {kon,koff}.

sol= ParametricNDSolve[eqox ,{x1, x2,x3},{t, 0, 100} ,{kon,koff}]

2.-Now it is computed  ∇(f(t), {a, b}) = {((df(t))/da),((df(t))/db)},

fa[a1_?NumberQ, b_?NumberQ, t_?NumberQ]:=D[x1[a,b],a][t]/.a->a1/.sol

fb[a_?NumberQ, b1_?NumberQ, t_?NumberQ]:=D[x2[a,b],b][t]/.b->b1/.sol

3.- Here is defined the number of points n to be used in the optimal design.

4.- It is evaluated ∇(f(t), β) at points t:{(Subscript[t, 0]),..., (Subscript[t, n])}, obtaining X = {(Subscript[X, 1]),..., (Subscript[X, p])} with (Subscript[X, 1]) = {((df((Subscript[t, 0])))/(Subscript[dβ, 1])), ...,((df((Subscript[t, n])))/(Subscript[dβ, 1])) }, ..., (Subscript[X, p]) = {((df((Subscript[t, 0])))/(Subscript[dβ, p])), ...,((df((Subscript[t, n])))/(Subscript[dβ, p])) }

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

X1[0.7, 0.11,  0.5] // AbsoluteTiming

Test OK( The same value that using Method 1)

6.- A typical election for compute the covariance matrix is assumed that  that the relationship between samples decays exponentially with increasing time-distance between them, that is  Γ = {(Subscript[l, ij])} with (Subscript[l, ij])= exp {ρ|(Subscript[t, j]) -(Subscript[t, j])|}.For computational purpose we have found more appropriate to use the distance (Subscript[d, i]) = (Subscript[t, i])-(Subscript[t, i-1]), instead of  (Subscript[t, i]), then (Subscript[t, i])= (∑_i)(Subscript[d, i]) being (Subscript[d, 0]) =(Subscript[t, 0]) . That is for a two points design . We suppose a 3-points design. The first is defined by the user 

 Γ where

FoldList[Plus, Subscript[d, 0], Table[(Subscript[d, i]),{i,n}]];
ff[i_,j_]:=Which[i==j,1,i<j,(E^(-ρ (∑_(k=i)^(j-1))(Subscript[d, k]))),i>j,(E^(-ρ (∑_(k=j)^(i-1))(Subscript[d, k])))];
nn=2;
Γ=Array[ff,{nn+1,nn+1}]

6.-  Now it is computed the convariance matrix Σ = (σ^2)  Γ

We take:

ρ=1;σ =1;

We will also  need give the initial values for the β parameters and the standard deviation of the measures. kon= 0.7, koff=0.11

7.- Then we can obtain the information matrix

M =  (X^T) (Σ^(-1)) X

m := X . Inverse[Σ].  Transpose[X];

m1[ti_] :=Transpose[Map[X1[0.7, 0.11, #]&,ti] ]. Inverse[Σ].Map[X1[0.7, 0.11, #]&,ti]

8.- Finally  the determinant of the information matrix is maximized as function of d0, d1 and d2. We constrain the d values to a maximun of (Subscript[d, i])=10 becouse to longer time the concentration will be very low (lower than the detection limit)

For n (number of observations)= 2

Γ={{1,(E^(-ρ d1))},{(E^(-ρ d1)),1}};

Σ = (σ^2)*Γ;

m1[ti_] :=Transpose[Map[X1[0.7, 0.11, #]&,ti] ]. Inverse[Σ].Map[X1[0.7, 0.11, #]&,ti]

obj[d0_?NumericQ,d1_?NumericQ]:=Det[m1[{d0,d1+d0}]]

sol2=NMaximize[ {obj[d0,d1],0.02<d0<10, 0.02<d1<10},{d0,d1}]

For n (number of observations)= 3

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2)))},{(E^(-ρ d1)),1,(E^(-ρ d2))},{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1}};

Σ = (σ^2)*Γ;

obj[d0_?NumericQ,d1_?NumericQ,d2_?NumericQ]:=Det[m1[{d0,d1+d0,d0+d1+d2}]]

sol3=NMaximize[ {obj[d0,d1,d2],0.02<d0<10,0.02<d1<10, 0.02<d2<10},{d0,d1,d2}]

For n (number of observations)= 4

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2))),(E^(-ρ (d1+d2+d3)))},
{(E^(-ρ d1)),1,(E^(-ρ d2)),(E^(-ρ (d2+d3)))},
{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1,(E^(-ρ d3))},
{(E^(-ρ (d1+d2+d3))),(E^(-ρ (d2+d3))),(E^(-ρ d3)),1}}; Σ = (σ^2)*Γ;

m1[ti_] :=Transpose[Map[X1[0.7, 0.11, #]&,ti] ]. Inverse[Σ].Map[X1[0.7, 0.11, #]&,ti]

obj[d0_?NumericQ,d1_?NumericQ,d2_?NumericQ, d3_?NumericQ]:=Det[m1[{d0,d1+d0,d0+d1+d2, d0+d1+d2+d3}]]

sol4=NMaximize[ {obj[d0,d1,d2, d3],0.02<d0<10,0.02<d1<10, 0.02<d2<10, 0.02<d3<10},{d0,d1,d2,d3}]

Conclusion: The observations will be taken: n, t0, t1,t2,t3}

{{"Observations 2:", d0,d1+d0}/. sol2[[2]],{"Observations 3:", d0,d1+d0,d0+d1+d2}/. sol3[[2]],{"Observations 4:", d0,d1+d0,d0+d1+d2, d0+d1+d2+d3}/. sol4[[2]]}

Quit[]

Michaelis-Menten

The model

In this case the drug transference between compartments will be considered as a linear kinetic process described by the transfer coefficients k12 and k21. However, the elimination process will be non-linear as it happens for instance in hepatic metabolism, and the elimination rate of the drug can mathematically be expressed by the Michaelis-Menten equation with parameters Vm=maximum transformation speed and km=Michaelis-Menten constant. The drug administration will be assumed to be an impulsive input (bolus).

1.- It is defined a model  f(t, β) where the unkown parameters are β = {a,b}.

Vmax = 0.2; km = 0.3

Vmax = 0.2; km = 0.3; V1=1;

eq1 = x1'[t]==-(Subscript[k, 12]) x1[t]+(Subscript[k, 21]) x2[t]-((Vmax x1[t]/V1)/(km+x1[t]/V1));

eq2 = x2'[t]==(Subscript[k, 12]) x1[t]-(Subscript[k, 21]) x2[t];

It can be solved as function of k12 and k21:

sol =ParametricNDSolve[{eq1, eq2,x1[0] == 1,x2[0] == 0},{x1,x2}, {t,0,100},{(Subscript[k, 12]),(Subscript[k, 21])}]

We wish  find  t :{(Subscript[t, 0]),..,  (Subscript[t, i]), ...  (Subscript[t, n])} )  using D-optimal design when the analitycal expresion of (Subscript[x, 1])(t,  (Subscript[k, a]),(Subscript[k, b])) can not be found. [f(t, β) =(Subscript[x, 1])(t,  (Subscript[k, a]),(Subscript[k, b]))]

2.-Now it is computed  ∇(f(t), {a, b}) = {((df(t))/da),((df(t))/db)},

fa[a1_?NumberQ, b_?NumberQ, t_?NumberQ]:=D[x1[a,b],a][t]/.a->a1/.sol

fb[a_?NumberQ, b1_?NumberQ, t_?NumberQ]:=D[x1[a,b],b][t]/.b->b1/.sol

3.-It is defined the number of points n to be used in the optimal design.

4.- It is evaluated ∇(f(t), β) at points t:{(Subscript[t, 0]),..., (Subscript[t, n])}, obtaining X = {(Subscript[X, 1]),..., (Subscript[X, p])} with (Subscript[X, 1]) = {((df((Subscript[t, 0])))/(Subscript[dβ, 1])), ...,((df((Subscript[t, n])))/(Subscript[dβ, 1])) }, ..., (Subscript[X, p]) = {((df((Subscript[t, 0])))/(Subscript[dβ, p])), ...,((df((Subscript[t, n])))/(Subscript[dβ, p])) }

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

We test using typical values of ki: k12 = 0.03, k21 = 0.02 for t:{1,3,10}

Map[X1[0.03, 0.02, #]&,{1,3,10}]//AbsoluteTiming

6.- A typical election for compute the covariance matrix is assumed that  that the relationship between samples decays exponentially with increasing time-distance between them, that is  Γ = {(Subscript[l, ij])} with (Subscript[l, ij])= exp {ρ|(Subscript[t, j]) -(Subscript[t, j])|}.For computational purpose we have found more appropriate to use the distance (Subscript[d, i]) = (Subscript[t, i])-(Subscript[t, i-1]), instead of  (Subscript[t, i]), then (Subscript[t, i])= (∑_i)(Subscript[d, i]) being (Subscript[d, 0]) =(Subscript[t, 0]) . That is for a two points design . We suppose a 3-points design. The first is defined by the user 

 Γ where

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2)))},{(E^(-ρ d1)),1,(E^(-ρ d2))},{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1}};;

6.-  Now it is computed the convariance matrix Σ = (σ^2)  Γ

Σ = (σ^2)*Γ;

We take:

ρ=1;σ =1;

Σ

We will also  need give the initial values for the β parameters and the standard deviation of the measures. ka= 0.03, kb=0.02

7.- Then we can obtain the information matrix

M =  (X^T) (Σ^(-1)) X

m := X . Inverse[Σ].  Transpose[X];

m1[ti_] :=Transpose[Map[X1[0.03, 0.02, #]&,ti] ]. Inverse[Σ].Map[X1[0.03, 0.02, #]&,ti]

obj[d0_?NumericQ,d1_?NumericQ,d2_?NumericQ]:=Det[m1[{d0,d1+d0,d0+d1+d2}]]

sol1=NMaximize[ {obj[d0,d1,d2],0<d0<10,0<d1<10, 0<d2<10},{d0,d1,d2}]//Timing

{d0,d1+d0,d0+d1+d2}/. sol1[[2,2]]

2.-Now it is computed  ∇(f(t), {a, b}) = {((df(t))/da),((df(t))/db)},

fa[a1_?NumberQ, b_?NumberQ, t_?NumberQ]:=D[x1[a,b],a][t]/.a->a1/.sol

fb[a_?NumberQ, b1_?NumberQ, t_?NumberQ]:=D[x1[a,b],b][t]/.b->b1/.sol

3.- Here is defined the number of points n to be used in the optimal design.

4.- It is evaluated ∇(f(t), β) at points t:{(Subscript[t, 0]),..., (Subscript[t, n])}, obtaining X = {(Subscript[X, 1]),..., (Subscript[X, p])} with (Subscript[X, 1]) = {((df((Subscript[t, 0])))/(Subscript[dβ, 1])), ...,((df((Subscript[t, n])))/(Subscript[dβ, 1])) }, ..., (Subscript[X, p]) = {((df((Subscript[t, 0])))/(Subscript[dβ, p])), ...,((df((Subscript[t, n])))/(Subscript[dβ, p])) }

X1[a_,b_, ti_]:={ fa[a,b, ti], fb[a,b, ti]}

We test using typical values of ki: ka = 0.03, kb = 0.02

Map[X1[0.03, 0.02, #]&,{1,3,10}]//AbsoluteTiming

6.- A typical election for compute the covariance matrix is assumed that  that the relationship between samples decays exponentially with increasing time-distance between them, that is  Γ = {(Subscript[l, ij])} with (Subscript[l, ij])= exp {ρ|(Subscript[t, j]) -(Subscript[t, j])|}.For computational purpose we have found more appropriate to use the distance (Subscript[d, i]) = (Subscript[t, i])-(Subscript[t, i-1]), instead of  (Subscript[t, i]), then (Subscript[t, i])= (∑_i)(Subscript[d, i]) being (Subscript[d, 0]) =(Subscript[t, 0]) . That is for a two points design . We suppose a 3-points design. The first is defined by the user 

 Γ where

Γ={{1,(E^(-ρ d1)),(E^(-ρ (d1+d2)))},{(E^(-ρ d1)),1,(E^(-ρ d2))},{(E^(-ρ (d1+d2))),(E^(-ρ d2)),1}};;

6.-  Now it is computed the convariance matrix Σ = (σ^2)  Γ

Σ = (σ^2)*Γ;

We take:

ρ=1;σ =1;

Σ

We will also  need give the initial values for the β parameters and the standard deviation of the measures. ka= 0.03, kb=0.02

7.- Then we can obtain the information matrix

M =  (X^T) (Σ^(-1)) X

m := X . Inverse[Σ].  Transpose[X];

m1[ti_] :=Transpose[Map[X1[0.03, 0.02, #]&,ti] ]. Inverse[Σ].Map[X1[0.03, 0.02, #]&,ti]

sol1=NMaximize[ {obj[d0,d1,d2],0<d0<10,0<d1<10, 0<d2<10},{d0,d1,d2}]

The observation should be taken (i):

{d0,d1+d0,d0+d1+d2}/. sol1[[2]]