Velten K. Mathematical Modeling and Simulation: Introduction for Scientists and Engineers

Подождите немного. Документ загружается.

3.10 More Examples 215

it is not really surprising that such a fundamental relation from enzyme kinetics

is used here since a number of enzyme-mediated reactions is involved in the

metabolization of nitrogen by the yeast cells [124].

Equation 3.284 is appropriate under largely isothermal conditions. Since we want

to consider the temperature as a (nonconstant) control variable here, the effects

of the temperature on μ

max

need to be taken into account, which can be done as

follows [122, 124, 125]:

max

(T) = 0.18 · exp



14 200 ·

T − 300

300RT



− 0.0054 · exp



121 000 ·

T − 300

300RT



(3.286)

Here, R = 8.314472 J · K

−1

· mol

−1

is the universal gas constant. The death

constant k

depends on the ethanol concentration E(t). In [121], the following

expression is used

= k · E (3.287)

The last equations describe the dynamics of the viable yeast cell concentration,

X(t). In these equations, we needed the nitrogen concentration, N(t), and the

ethanol concentration, E(t). This means we need equations characterizing N(t)and

E(t). Let us begin with an ODE describing the dynamics of E(t).

3.10.2.3 Ethanol and Sugar

Cramer describes the ethanol production rate proportional to the available amount

of yeast cells, that is,

= βX (3.288)

where, similar to Equation 3.284, the speciﬁc ethanol production rate β depends

on the available sugar concentration, S(t):

β = β

max

+ S

(3.289)

Again, the last two equations hold for essentially isothermal conditions, and a

temperature dependence of β

max

needs to be taken into account here. Following

[122], we use an expression derived from data in [121, 126]:

max

(T − 273.15) = β

max,24

◦

· (0.00132 ·T

+ 0.00987 · T − 0.00781)

(3.290)

It remains to describe N(t)andS(t) that are used in the last equations. As discussed

above, the dynamics of the sugar concentration S(t ) is intimately related with the

216 3 Mechanistic Models I: ODEs

ethanol production by the yeast cells. This is reﬂected by the fact that an ODE very

similar to Equation 3.288 can be used to describe the sugar dynamics [121]:

=−

X (3.291)

In this equation, the coefﬁcient Y

describes the stoichiometric bioconversion

of sugar into ethanol.

3.10.2.4 Nitrogen

Remember that we wanted to use N(t ) as one of our control variables. In practice,

any kind of ﬁne tuning of N(t) is hard to achieve. Even the measurement of nitrogen

levels – which would be an essential prerequisite of controlling N(t) levels – is a

very difﬁcult task [122, 124]. Nevertheless, nitrogen is an important variable in

this process and people try to make sure that there is sufﬁcient nitrogen in the

fermenter. Since they usually do not know the actual nitrogen level N(t), this

is done largely in a ‘‘black box fashion’’, which means that nitrogen is added

to the fermenter at times which have proven to be reasonable based on the

experience made in previous fermentations. Let us denote with N

add

the nitrogen

concentrations added at times t

add

(i = 1, ..., n). From the process engineers point

of view, these are the nitrogen control variables, and they give rise to an essentially

unknown dynamics of N(t). One of the beneﬁts of the model described here is that

based on the nitrogen levels computed by the model, the process engineer can at

least gain a qualitative idea of the effects that the nitrogen additions have on N(t).

Note 3.10.1 (Qualitative optimization) In a situation where an exact quan-

titative determination of a state variable cannot be achieved, mechanistic

mathematical models often provide information about its qualitative behav-

ior, which can be used for optimization. Such a ‘‘qualitative optimization’’ of

a process is a deﬁnite step beyond pure ‘‘black box optimization’’ based on

phenomenological models.

Since N(t) is required in the above equations, it must be computed based on the

nitrogen additions N

add

at times t

add

. To do this, we start with Cramer’s nitrogen

balance equation as follows:

=−

· X (3.292)

This equation is based on the assumption that nitrogen is consumed proportional

to the growth rate of yeast. The yield coefﬁcient Y

is used to convert from nitrogen

to yeast biomass. Now the nitrogen additions N

add

at times t

add

must be incorporated

into the equations. Unfortunately, R’s

lsoda does not provide us with an option that

would allow us to impose a discontinuous change of N(t)byanamountofN

add

times t

add

. Principally, this could be done by splitting the overall integration interval

3.10 More Examples 217

[0, T]intosubintervals[0,t

add

], [t

add

, t

add

], and so on, applying lsoda separately to

each of the subintervals and using appropriate settings of the initial values that

account for the nitrogen additions N

add

. In our case, a more elegant and sufﬁciently

accurate alternative is the assumption that the nitrogen additions are not given

instantaneously, but at a constant rate over 1 h. This turns Equation 3.292 into

= r(t) −

· X (3.293)

where

r(t) =



add

if t

add

−

< t < t

add

, i = 1, ..., n

0 else

(3.294)

The overall model is a system of four ODEs that can be summarized as follows:

Wine fermentation model

= μ(T, N) · X − k · E · X (3.295)

= r(t) −

· μ(T, N) ·X (3.296)

= β(T, S) · X (3.297)

=−

β(T, S) ·X (3.298)

We have used the notations μ(T, N)andβ(T, S) to emphasize the nonlinear

dependence of these parameters on their arguments as discussed above.

3.10.2.5 Using a Hand-ﬁt to Estimate N

Now let us see if the above model is capable to reproduce the concrete fermentation

data in Figure 3.18. In this fermentation, the temperature has been varied between

12 and 18

◦

C by the process engineers (Figure 3.18), and nitrogen has been added

at two times during the fermentation (see t

add

and N

add

in Table 3.4). As a result of

these measures, the ﬁgure shows data of X(t), E (t)andS(t). As it was mentioned

above, the measurement of N(t) is difﬁcult and this is why no nitrogen data are

available from this trial. Figure 3.18 was produced using the program Fermentation.r

in the book software (set

PlotData=TRUE in the ‘‘Program Control Parameters’’

section).

Fermentation.r will now also be used to solve the ODE system, Equations

3.295–3.298, and to estimate parameters of these equations from the data. This

code has been derived by an appropriate editing of

ODEFitEx1.r (Section 3.9.1).

Some data handling aspects of

Fermentation.r closely related to the parameter

218 3 Mechanistic Models I: ODEs

0 100

t (h)

200

100

150

Temperature (°C)

Sugar (g l

−1

)

3.5

2.5

100

Cell biomass (g l

−1

)

Ethanol (g l

−1

)

1.5

0.5

0 50 100 150 200 250 300 350 0 100 200 300

200 300

0 100

t (h)

t (h) t (h)

200 300

Fig. 3.18 Wine fermentation data from [122]. See fermentation.csv in the book software.

estimation procedure have already been discussed in Section 3.9.2.2. Other aspects

of that code that are more closely related with the wine fermentation model

are addressed below. See also the discussion of

ODEEx1.r and ODEEx2.r in

Section 3.8.3 for more information on the general way how you should work with

this kind of codes.

Fermentation.r is based on Equations 3.295–3.298 and uses

the parameter values given in Table 3.4.

Looking at the table, you will note that there are two unknown parameters: k

and N

. N

is hard to measure because it refers to the yeast available nitrogen

that cannot be easily characterized experimentally [122, 128]. Regarding k, [121]

suggests a value of k = 10

−4

l/g ethanol/h for a temperature of 24

◦

C. This value

could have been used here as a ﬁxed value in the simulations, but [122] says it

is better to estimate k from the data since this parameter is strongly temperature

dependent, and the temperatures in the trial investigated here are substantially

below 24

◦

C (Figure 3.18). Therefore, both N

and k are estimated from the data in

Fermentation.r using the procedure explained in Section 3.9.

First of all, starting values of k and N

are needed as initial values for the

parameter estimation procedure, that is, values in the right order of magnitude.

As explained above, one can try to get such approximate values from the literature

and/or by a hand-ﬁtting of the parameters until the solutions of the ODE system

are at least in the neighborhood of the data. Regarding k, we can use the literature

value k = 10

−4

mentioned above. Since we do not have any information on N

the hand-ﬁtting procedure is used for this parameter. A hand-ﬁt of N

(setting

NonlinRegress=FALSE and PlotStateData=TRUE in the ‘‘program control pa-

rameters’’ Section of

Fermentation.r) leads us, for example, to Figure 3.19,

which is obtained for k = 10

−4

and N

= 0.1.

3.10 More Examples 219

0 100

t (h)

200

100

150

3.5

2.5

1.5

0.5

Sugar (g l

−1

)

Cell biomass (g l

−1

)

100

0.08

0.04

0.00

Ethanol (g l

−1

)

Nitrogen (g l

−1

)

0 100 200 300 400 0 100 200 300 400

200 300 400 0 100

t (h)

t (h) t (h)

200 300 400

= 0.751

= 0.743

= 0.566

Fig. 3.19 Solution of Equations 3.295–3.298 using k = 10

−4

= 0.1 and the parameters in Table 3.4.

3.10.2.6 Parameter Estimation

The ﬁne tuning of the parameters to be estimated can then be performed based on

R’s

nls function as discussed in Section 3.9, using k = 10

−4

and the hand-ﬁtted

= 0.1asstartingvaluesofk and N

. Setting NonlinRegress=TRUE in the

‘‘program control parameters’’ section of

Fermentation.r, the result shown in

Figure 3.20 is obtained. As can be seen, a very good ﬁt between the solution of

the ODE system and the data is obtained, with all R

values above 0.97. This is

a very good result since this is a ﬁt between three nonlinear curves and data that

was obtained by a tuning of only two parameters. If you have achieved a result

like this, the next step is to test the predictive power of your model, applying it

to data not used for the parameter estimation. This has been done by Blank in

[122]. Blank shows that this model performs very well on unknown data in several

cases, although there are also datasets where the results are unsatisfactory – which

is not really surprising considering the aforementioned complexity of the wine

fermentation process. From the point of view of mathematical modeling, the latter

datasets producing the unsatisfactory results are the interesting ones, because they

give us hints for further improvements of the model.

Note that the relative and absolute tolerances

rtolDef and atolDef have both

been set to 10

−5

in Fermentation.r. This is based on the ‘‘rule of thumb’’

explained in Section 3.8.3.5: Considering the fact that the computations involve

nitrogen additions of 0.03 gl

−1

(see t

add

and N

add

in Table 3.4), it is clear that

atolDef must be smaller than 10

−2

. For example, we can set it to 10

−3

.Regarding

rtolDef, Figure 3.18 shows that the sugar and ethanol data have at least three

220 3 Mechanistic Models I: ODEs

0 100

t (h)

200

100

150

3.5

2.5

1.5

0.5

Sugar (g l

−1

)

Cell biomass (g l

−1

)

100

0.10

0.15

0.05

0.00

Ethanol (g l

−1

)

Nitrogen (g l

−1

)

0 100 200 300 400 0 100 200 300 400

200 300 400 0 100

t (h)

t (h) t (h)

200 300 400

= 0.996

= 0.984

= 0.976

Fig. 3.20 Solution of Equations 3.295–3.298 using the

estimated coefﬁcients k = 6.637202 × 10

−5

and N

1.678096 × 10

−1

, and the parameters in Table 3.4.

signiﬁcant digits, which means we should set rtolDef also, for example, to 10

−3

Fermentation.r, both tolerances were set to 10

−5

,thatis,twoordersof

magnitude smaller, which was based on the heuristical procedure for the choice

rtolDef and atolDef explained in Section 3.8.3.5 (a further decrease in the

tolerances below 10

−5

did not affect the result anymore).

3.10.2.7 Problems with Nonautonomous Models

Finally, let us look through Fermentation.r,focusingonnewaspectsinthis

code. As it was said above,

Fermentation.r was derived from ODEFitEx1.r

which was discussed in Section 3.9.1. Within Fermentation.r,theODEsystem

Equations 3.295–3.298 is deﬁned in a function

dfn as follows:

1: dfn <-

2: function(t, y, p)

3: {

4: X=y[1]

5: N=y[2]

6: E=y[3]

7: S=y[4]

8: k=p[1]

9: dXdt=mu(T(t),N)*X-k*E*X

10: dNdt= r(t)-mu(T(t),N)*X/Yxn

11: dEdt=beta(T(t),S)*X

(3.299)

3.10 More Examples 221

12: dSdt=-beta(T(t),S)*X/Yes

13: list(c(dXdt,dNdt,dEdt,dSdt))

14: }

The general structure of the dfn function has been discussed in Section 3.8.3.

The correspondence between Equations 3.295–3.298 with lines 9–12 of program

3.299 is obvious. You may be irritated by the fact that the time dependence of

r(t) and T(t) is written explicitly in the code, while this is not done for the state

variables

X, N, E and S. This is related with the fact that Equations 3.295–3.298 are

a nonautonomous system of ODEs. As it was discussed in Section 3.5.3, this means

that at least one of the right-hand sides of these equations depends on t not only

(implicitly) via the state variables, but also via some explicitly given functions of t.

Looking at the equations, you see that all right-hand sides of the system depend

explicitly on the temperature T(t) or on the nitrogen supply rate r(t). When writing

down the

dfn function for nonautonomous ODEs, you must take care that only the

state variables are written without an explicit ‘‘

(t)’’, while all other time-dependent

function must be written similar to

T(t), r(t) and so on.

But even if you have done this correctly, you might run into another problem

characteristic of nonautonomous ODEs when you try to compute the numerical

solution. Remember our deﬁnition of the nitrogen supply rate in Equation 3.294.

In the above example, nitrogen is supplied at two times only, t

add

= 130 h and

add

= 181 h. Since the overall time interval in the example is from 0 to 400 h,

this means that the function r(t) deﬁned in Equation 3.294 is almost everywhere

zero, except for two small 1-h intervals around t

add

= 130 h and t

add

= 181 h.

Remember also the discussion of

lsoda in Section 3.8: lsoda determines its

stepsize h automatically, choosing the stepsize as small as it is required by the error

tolerances, but not smaller. Now if

lsoda sets the stepsize relatively large around

add

= 130 h and t

add

= 181 h, it may happen that it ignores the nitrogen supply

rate r(t) in the sense that r(t) is never evaluated at one of its nonzero points. Then,

although you prescribe a substantial nitrogen supply in two small 1-h intervals

around t

add

= 130 h and t

add

= 181 h, you will not see this in the solution, that is,

the problem is solved by

lsoda in a way as if there would be no nitrogen supply

at all (r(t) = 0). Implementing the fermentation model yourself by an appropriate

editing of

ODEFitEx1.r, you will see that you get exactly this problem when you

do not apply appropriate measures to avoid this.

Note 3.10.2 (Problem with short external events) Numerical solutions of

nonautonomous ODE systems may ignore short external events (such as nitrogen

addition in the wine fermentation model).

The ﬁrst and simplest thing that one can do is to prescribe a maximum stepsize

lsoda, which can be done using an argument of lsoda called hmax (see [108,

109] and R’s hep pages). Choosing

hmax sufﬁciently small, one can make sure that

r(t) is used in the computation as desired. For example,, based on

hmax=1/10, r(t)

222 3 Mechanistic Models I: ODEs

will be evaluated 10 times during each of those 2 h of nitrogen supply. Similar to

the heuristics for an optimal choice of the stepsize described above (Note 3.8.1),

one can then optimize the size of

hmax, making it small enough such that a further

decrease of

hmax does not affect your results, but not smaller since this would

further increase your computation times, roundoff errors and so on as discussed

in Section 3.8.

A second option that can avoid this unwanted ‘‘blindness’’ of

lsoda with respect

to r(t) lies in an appropriate choice of the time vector supplied in the call of

lsoda.

As it was mentioned in the above discussion of

ODEEx1.r,thetimevectorin

the call of

lsoda (line 2 in program 3.249) deﬁnes the times at which lsoda is

expected to generate numerical approximations of the state variables. Now if we

want

lsoda to ‘‘see’’ those 2 h where the nitrogen supply rate r(t) is different from

zero, we can deﬁne that time vector such that it involves several times within those

2 h. Exactly this has been done in

Fermentation.r, using a time vector supplied

lsoda called tInt. Check the deﬁnition of tInt in Fermentation.r to see

that it really does what has been described above.

3.10.2.8 Converting Data into a Function

The temperature T(t) in Equations 3.295–3.298 needs some special treatment since

it is not known as a mathematical function, but rather in the form of experimental

data (columns 1–2 in

fermentation.csv in the book software). The conversion

of such experimental data into a form that can be treated by software such as

lsoda

is a standard problem in the numerical treatment of ODEs, and it is usually solved

by a simple linear interpolation of the data. Denoting with (t

, T

)(i = 1, ..., n)the

experimental data, this means that T(t) is deﬁned in each of the intervals [t

i−1

, t

]

(i = 2, ..., n) as follows:

T(t) = T

i−1

t − t

i−1

− t

i−1

· (T

− T

i−1

) (3.300)

3.10.2.9 Using Weighting Factors

The wine fermentation model is an example where it makes sense to use the

weighting factors w

in the residual sum of squares, Equation 3.260 (Section 3.9.2).

As it was discussed there, weighting factors should be used in situations where

the measurement data cover several orders of magnitude. As Figure 3.18 shows,

the sugar and cell biomass data differ by about two orders of magnitude, and the

sugar and ethanol data cover two orders of magnitude. In [41] it is suggested that

aweightingw

= 1/σ

j = 1, ..., m

should be used in this situation, where σ

the standard deviation of the measurement data referring to state variable i. [122]

reports the following values for the state variables of the wine fermentation model:

= 0.1, σ

= 0.7, σ

= 0.7. This gives the weighting factors

= 100 (3.301)

3.10 More Examples 223

≈ 2 (3.302)

≈ 2 (3.303)

Here, w

refers to all w

in the residual sum of squares, Equation 3.260, where

i refers to X (similar for w

and w

). In Fermentation.r, this weighting is

prescribed in the ‘‘program control parameters’’ section as follows:

WhichData=c(1,1,50)

(3.304)

which means that X (third entry in WhichData) receives 50 times more weighting

compared to S and E (ﬁrst and second entries in

WhichData). Looking into the call

nls further below in Fermentation.r,youcanseehowWhichData is used to

deﬁne an option of

nls called weights.

3.10.3

Pharmacokinetics

Most therapeutic drugs are ingested orally and then enter the blood stream via

gastrointestinal (GI) absorption. For a therapeutic drug to be effective, it is of

course important to know how fast the drug enters the blood stream, and how the

drug concentration changes depending on time. In (probably) the simplest possible

approach, we would just consider two state variables: G(t), the concentration of the

drug in the GI tract, and B(t), the concentration of the drug in the blood (we assume

units of micrograms per milliliter). If the application rate D(t) (micrograms per

milliliter per hour) expresses the drug dosage regime as seen by the GI tract, the

dynamics of the system can be described by the following ODE system [116, 129]:



(t) =−aG(t) +D (t) (3.305)



(t) = aG(t) − bB(t) (3.306)

Equation 3.306 says that the blood concentration B(t) grows proportionally to the

concentration in the GI tract, G(t), and it decays proportionally to B(t), where a

and b (per hour) are the constants of proportionality. The term aG(t) that describes

the increase of the blood concentration of the drug by absorption appears with an

inverse sign in Equation 3.305, expressing the mass balance of the drug in the

sense that precisely the amount of drug that enters the blood by absorption is lost

in the GI tract. Of course, the dosage D(t) causes an increase of G(t)andthus

appears with a ‘‘+’’-sign in Equation 3.305.

To show the qualitative behavior of the solution, we use the parameter settings

suggested in [116]: G(0) = B(0) = 0 (i.e. no drug has been applied before t = 0),

a = ln(2)/2, b = ln(2)/5, and

D(t) = 2 ·



n=0



H (t − 6n) − H



t −



6n +



(3.307)

224 3 Mechanistic Models I: ODEs

where H(x) is the Heaviside step function:

H (x) =



0ifx < 0

1ifx ≥ 0

(3.308)

Equation 3.307 yields the drug dosage regime shown in Figure 3.21, which can

be thought of as representing a situation where the drug is applied every 6 h, for

example, at 6 a.m., 12 a.m., and 6 p.m. As can be seen in the ﬁgure (and as

it follows from Equation 3.307), it is assumed here that the drug is supplied at

a continuous rate of 2 μgml

−1

, and that this rate is held constant for 1/2 h

after drug ingestion. Figure 3.22 shows the resulting patterns of B(t)andG(t)that

are obtained by a solution of the model, Equations 3.305 and 3.306. To solve the

mathematical problem, the code

ODEEx3.r hasbeenused.Thiscodeisapartof

the book software (see Appendix A), and it has been obtained by an editing of

ODEEx2.r, a similar code that has been discussed in Section 3.8.3 above.

t (h)

0 5 10 15 20 25 3

D (μg ml

−1

)

0.0

0.5

1.0

1.5

2.0

Fig. 3.21 Assumed pattern of drug dosage D(t) corresponding to Equation 3.307.

0510

0.0

0.2

0.4

0.6

15 20 25 30

t (h)

G (μg ml

−1

)

0510

0.0

0.4

0.8

1.2

15 20 25 30

t (h)

(a) (b)

B (μg ml

−1

)

Fig. 3.22 Drug concentration (a) in the gastrointestinal tract and (b) in the blood.