104,99 €
Mehr erfahren.
- Herausgeber: John Wiley & Sons
- Kategorie: Wissenschaft und neue Technologien
- Sprache: Englisch
Features a solid foundation of mathematical and computational tools to formulate and solve real-world ODE problems across various fields With a step-by-step approach to solving ordinary differential equations (ODEs), Differential Equation Analysis in Biomedical Science and Engineering: Ordinary Differential Equation Applications with R successfully applies computational techniques for solving real-world ODE problems that are found in a variety of fields, including chemistry, physics, biology, and physiology. The book provides readers with the necessary knowledge to reproduce and extend the computed numerical solutions and is a valuable resource for dealing with a broad class of linear and nonlinear ordinary differential equations. The author's primary focus is on models expressed as systems of ODEs, which generally result by neglecting spatial effects so that the ODE dependent variables are uniform in space. Therefore, time is the independent variable in most applications of ODE systems. As such, the book emphasizes details of the numerical algorithms and how the solutions were computed. Featuring computer-based mathematical models for solving real-world problems in the biological and biomedical sciences and engineering, the book also includes: * R routines to facilitate the immediate use of computation for solving differential equation problems without having to first learn the basic concepts of numerical analysis and programming for ODEs * Models as systems of ODEs with explanations of the associated chemistry, physics, biology, and physiology as well as the algebraic equations used to calculate intermediate variables * Numerical solutions of the presented model equations with a discussion of the important features of the solutions * Aspects of general ODE computation through various biomolecular science and engineering applications Differential Equation Analysis in Biomedical Science and Engineering: Ordinary Differential Equation Applications with R is an excellent reference for researchers, scientists, clinicians, medical researchers, engineers, statisticians, epidemiologists, and pharmacokineticists who are interested in both clinical applications and interpretation of experimental data with mathematical models in order to efficiently solve the associated differential equations. The book is also useful as a textbook for graduate-level courses in mathematics, biomedical science and engineering, biology, biophysics, biochemistry, medicine, and engineering.
Sie lesen das E-Book in den Legimi-Apps auf:
Seitenzahl: 360
Veröffentlichungsjahr: 2014
Ähnliche
Table of Contents
Cover
Title Page
Copyright
Dedication
Preface
Chapter 1: Introduction to Ordinary Differential Equation Analysis: Bioreactor Dynamics
1.1 Introduction
1.2 A ODE System for a Bioreactor
1.3 In-Line ODE Routine
1.4 Numerical and Graphical Outputs
1.5 Separate ODE Routine
1.6 Alternative Forms of ODE Coding
1.7 ODE Integrator Selection
1.8 Euler Method
1.9 Accuracy and Stability Constraints
1.10 Modified Euler Method as a Runge–Kutta Method
1.11 Modified Euler Method as an Embedded Method
1.12 Classic Fourth-Order Runge–Kutta Method as an Embedded Method
1.13 RKF45 Method
References
Chapter 2: Diabetes Glucose Tolerance Test
2.1 Introduction
2.2 Mathematical Model
2.3 Computer Analysis of the Mathematical Model
2.4 Conclusions
References
Chapter 3: Apoptosis
3.1 Introduction
3.2 Mathematical Model
3.6 Base Case with Variation in ICs
3.7 Variation in ODEs
3.8 Selection of Units
3.9 Model Solution with RKF45
3.10 Conclusion
Reference
Chapter 4: Dynamic Neuron Model
4.1 Introduction
4.2 The Dynamic Neuron Model
4.3 ODE Numerical Integration
4.4 Conclusions
References
Chapter 5: Stem Cell Differentiation
5.1 Introduction
5.2 Model Equations
5.3 R Routines
5.4 Summary
Reference
Chapter 6: Acetylcholine Neurocycle1
6.1 Introduction
6.2 ODE Model
6.3 Numerical Solution of the Model
6.4 Model Output
6.5 ODE/PDE Model
Appendix A1: IC Vector by a Differential Levenberg Marquardt Method
References
Chapter 7: Tuberculosis with Differential Infectivity1
7.1 Introduction
7.2 Mathematical Model
7.3 R Routines for the ODE Model
7.4 Model Output
7.5 Conclusions
References
Chapter 8: Corneal Curvature
8.1 Introduction
8.2 Model Equations
8.3 Method of Lines Solution
8.4 R Routines
8.5 Numerical Solution
8.6 Error Analysis of the Numerical Solution
8.7 Library Routines for Differentiation in Space
8.8 Summary
References
Appendix A1: Stiff ODE Integration
A1.1 Introduction
A1.2 Analytical Solution of Second-Order ODE System
A1.3 Eigenvalue Stability Analysis
A1.4 BDF Methods for Stiff ODEs
A1.5 R Program for First-Order BDF Method
A1.6 Numerical Output from the BDF Integration
A1.7 Alternative Programming of the BDF Integration
A1.8 Second-Order BDF Integration
A1.9 Third-Order BDF Integration
A1.10 Conclusions
References
Index
Pages
ix
x
xi
xii
xiii
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
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
125
126
127
128
129
130
124
131
132
133
134
135
136
137
138
139
140
141
142
143
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
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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
236
237
238
239
241
242
243
245
244
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
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
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
417
418
419
420
421
422
423
Guide
Cover
Table of Contents
Preface
Chapter 1: Introduction to Ordinary Differential Equation Analysis: Bioreactor Dynamics
List of Illustrations
Figure 1.1
Figure 1.2
Figure 2.1
Figure 2.2
Figure 2.3
Figure 2.4
Figure 3.1
Figure 3.2
Figure 3.3
Figure 3.4
Figure 4.1
Figure 4.3
Figure 4.2
Figure 4.4
Figure 4.5
Figure 5.1
Figure 5.2
Figure 5.3
Figure 5.4
Figure 6.1
Figure 6.2
Figure 6.6
Figure 6.3
Figure 6.7
Figure 6.10
Figure 6.11
Figure 7.1
Figure 7.2
Figure 7.3
Figure 8.1
Figure 8.3
Figure 8.2
Figure 8.4
Figure 8.6
List of Tables
Table 1.1
Table 1.2
Table 1.3
Table 1.4
Table 1.5
Table 1.6
Table 1.7
Table 1.8
Table 2.1
Table 2.2
Table 2.3
Table 2.4
Table 2.5
Table 3.1
Table 3.2
Table 3.3
Table 3.4
Table 3.5
Table 3.6
Table 3.7
Table 3.8
Table 4.1
Table 4.2
Table 4.3a
Table 4.3b
Table 5.1
Table 6.1
Table 6.2
Table 6.3
Table 6.4
Table 6.5a
Table 6.5b
Table 6.6a
Table 6.6b
Table 6.7a
Table 6.7b
Table 7.1
Table 7.2
Table 7.3
Table 7.4
Table 8.1
Table 8.2
Table 8.3
Table 8.4
Table A1.1
Table A1.2
Table A1.3
Table A1.4
Table A1.5
Table A1.6
Differential Equation Analysis in Biomedical Science and Engineering
Ordinary Differential Equation Applications with R
William E. Schiesser
Cover Design: Wiley
Cover Illustration: Courtesy of the author
Copyright © 2014 by John Wiley & Sons, Inc. All rights reserved.
Published by John Wiley & Sons, Inc., Hoboken, New Jersey.
Published simultaneously in Canada.
No part of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, scanning, or otherwise, except as permitted under Section 107 or 108 of the 1976 United States Copyright Act, without either the prior written permission of the Publisher, or authorization through payment of the appropriate per-copy fee to the Copyright Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923, (978) 750-8400, fax (978) 750-4470, or on the web at www.copyright.com. Requests to the Publisher for permission should be addressed to the Permissions Department, John Wiley & Sons, Inc., 111 River Street, Hoboken, NJ 07030, (201) 748-6011, fax (201) 748-6008, or online at http://www.wiley.com/go/permission.
Limit of Liability/Disclaimer of Warranty: While the publisher and author have used their best efforts in preparing this book, they make no representations or warranties with respect to the accuracy or completeness of the contents of this book and specifically disclaim any implied warranties of merchantability or fitness for a particular purpose. No warranty may be created or extended by sales representatives or written sales materials. The advice and strategies contained herein may not be suitable for your situation. You should consult with a professional where appropriate. Neither the publisher nor author shall be liable for any loss of profit or any other commercial damages, including but not limited to special, incidental, consequential, or other damages.
For general information on our other products and services or for technical support, please contact our Customer Care Department within the United States at (800) 762-2974, outside the United States at (317) 572-3993 or fax (317) 572-4002.
Wiley also publishes its books in a variety of electronic formats. Some content that appears in print may not be available in electronic formats. For more information about Wiley products, visit our web site at www.wiley.com.
Library of Congress Cataloging-in-Publication Data:
Schiesser, W. E.
Differential equation analysis in biomedical science and engineering : ordinary differential equation
applications with R / William E. Schiesser, Department of Chemical Engineering, Lehigh University,
Bethlehem, PA.
pages cm
Includes bibliographical references and index.
ISBN 978-1-118-70548-3 (cloth)
1. Biomedical engineering–Mathematics. 2. Biomathematics. 3. Bioreactors–Data processing. 4. Differential equations. I. Title.
R857.M34O34 2013
610.280285–dc23
2013020440
To John Backus and Leon Lapidus,
and their research groups
Preface
This book focuses on the rapidly expanding development and use of computer-based mathematical models in the life sciences, designated here as biomedical science and engineering (BMSE). The mathematical models are stated as systems of ordinary differential equations (ODEs) and generally come from papers in the current research literature that typically include the following steps:
The model is presented as a system of ODEs that explain associated chemistry, physics, biology, and physiology.
A numerical solution to the model equations is presented, particularly a discussion of the important features of the solution.
What is missing in this two-step approach are the details of how the solution was computed, particularly the details of the numerical algorithms. Also, because of the limited length of a research paper, the computer code used to produce the numerical solution is not provided. Thus, for the reader to reproduce (confirm) the solution and extend it is virtually impossible with reasonable effort.
The intent of this book is to fill in the steps for selected example applications that will give the reader the knowledge to reproduce and possibly extend the numerical solutions with reasonable effort. Specifically, the numerical algorithms are discussed in some detail, with additional background references, so that the reader will have some understanding of how the calculations were performed, and a set of transportable routines in R1 is provided so that the reader can execute to produce and extend the solutions.
Thus, the typical format of a chapter includes the following steps:
The model is presented as a system of ODEs with explanation of the associated chemistry, physics, biology, and physiology. The requirements of a well-posed set of equations such as the number of dependent variables, the number of ODEs, algebraic equations used to calculate intermediate variables, and the initial conditions for the ODEs are included (which is often not the case in research papers so that all of the details of the model are not included or known to the reader).
The features of the model that determine the selection of numerical algorithms are discussed; for example, how derivatives are approximated, whether the ODEs are nonstiff or stiff, and therefore, whether an explicit or implicit integration algorithm is used. The computational requirements of the particular selected algorithms are identified such as the solution of nonlinear equations, banded matrix processing, or sparse matrix processing.
The routines that are the programming of the ODEs and numerical algorithms are completely listed and then each section of code is explained, including referral to the mathematical model and the algorithms. Thus, all of the computational details for producing a numerical solution are in one place. Referring to another source for the software, possibly with little or no documentation, is thereby avoided.
A numerical solution to the model equations is presented, particularly a discussion of the important features of the solution.
The accuracy of the computed solution is inferred using established methods such as and refinement. Alternative algorithms and computational details are considered, particularly to extend the model and the numerical solution.
In this way, a complete picture of the model and its computer implementation is provided without having to try to fill in the details of the numerical analysis, algorithms, and computer programming (often a time-consuming procedure that leads to an incomplete and unsatisfactory result). The presentation is not heavily mathematical, for example, no theorems and proofs, but rather is in terms of detailed examples of BMSE applications.
End of the chapter problems have not been provided. Rather, the instructor can readily construct problems and assignments that will be in accordance with the interests and objectives of the instructor. This can be done in several ways by developing variations and extensions of the applications discussed in the chapters. For example,
Parameters in the model equations can be varied, and the effects on the computed solutions are observed and explained. Exploratory questions can be posed such as whether the changes in the solutions are as expected. In addition, the terms in the right- hand sides (RHSs) of the ODEs (without the derivatives in the independent variable) can be computed and displayed numerically and graphically to explain in detail why the parameter changes had the observed effect. The computation and display of ODE RHS terms is illustrated in selected chapters to serve as a guide.
Additional terms can be added to the ODE RHSs to model physical, chemical, and biological effects that might be significant in determining the characteristics of the problem system. These additional terms can be computed and displayed along with the original terms to observe which terms have a significant effect on the computed model solution.
One or more ODEs can be added to an existing model to include additional phenomena that are considered possibly relevant to the analysis and understanding of the problem system.
An entirely new model can be proposed and programmed for comparison with an existing model. The existing routines might serve as a starting point, for example, a template.
These suggested problem formats are in the order of increasing generality to encourage the reader to explore new directions, including the revision of an existing model and the creation of a new model. This process is facilitated through the availability of existing routines for a model that can first be executed and then modified. The trial-and-error development of a model can be explored, particularly if experimental data that can be used as the basis for model development are provided, starting from parameter estimation based on the comparison of experimentally measured data and computed solutions from an existing model up to the development of a new model to interpret the data.
The focus of this book is primarily on models expressed as systems of ODEs that generally result by neglecting spatial effects so that the dependent variables of ODE, for example, concentrations, are uniform in space. The independent variable of the ODE systems is, therefore, time in most applications.
The assumption of spatial uniformity is quite accurate for a spectrum of BMSE systems, for example, a mixed reactor as considered in Chapter 1. However, the assumption of spatial uniformity is not correct for BMSE systems that function according to the spatial distribution of their principal variables, such as corneal curvature discussed in the concluding Chapter 8. This chapter, therefore, serves as anintroduction to models based on partial differential equations (PDEs) in which both space and time are the independent variables. The additional spatial variables require boundary conditions for a complete specification of the PDE model. These details are introduced in Chapter 8 and are considered in detail in a companion volume titled Differential Equation Analysis in Biomedical Science and Engineering: Partial Differential Equation Applications with R.
In summary, my intention is to provide a set of basic computational procedures for ODE/PDE models in the two volumes that readers can use with modest effort without becoming deeply involved in the details of the numerical methods for ODE/PDEs and computer programming. All of the R routines discussed in the two volumes are available from a software download site, booksupport.wiley.com, which requires the ISBN: 9781118705483 for the ODE volume or 9781118705186 for the PDE volume. I welcome comments and will be pleased to respond to questions to the extent possible by e-mail ([email protected]).
William E. Schiesser
Bethlehem, PA
January 2014
1 R is a quality open source scientific programming system that can be easily downloaded from the Internet (http://www.R-project.org/). In particular, R has (i) vector-matrix operations that facilitate the programming of linear algebra, (ii) a library of quality ODE integrators, and (iii) graphical utilities for the presentation of the numerical solutions. All of these features and utilities are demonstrated through the applications in this book.
Chapter 1Introduction to Ordinary Differential Equation Analysis: Bioreactor Dynamics
1.1 Introduction
Mathematical models formulated as systems of ordinary differential equations (ODEs) and partial differential equations (PDEs) have been reported for a spectrum of applications in biomedical science and engineering (BMSE). The intent of this research is to provide a quantitative understanding of the biological, chemical, and physical phenomena that determine the characteristics of BMSE systems and to provide a framework for the analysis and interpretation of experimental data observed in the study of BMSE systems.
In the subsequent discussion in this chapter, we consider the programming of a (seven equations in seven unknowns) ODE system to illustrate the integration (solution) of ODE systems using R, a quality, open source, scientific programming system [10]. The intent is to provide the reader with a complete and thoroughly documented example of the numerical integration of an ODE system, including (i) the use of library ODE integrators, (ii) the programming of ODE integration algorithms, and (iii) graphical output of the numerical solutions. This example application can then serve as a prototype or template which the reader can modify and extend for an ODE model of interest.
1.2 A ODE System for a Bioreactor
The reaction system for the conversion of xylose to ethanol by fermentation is now formulated and coded (programmed) in R. The ODE model is discussed in detail in [[1], pp 35–42]; this discussion is recommended as a starting point for the details of the chemical reactions, particularly the various intermediates, so that the discussion to follow can concentrate on the numerical algorithms and R programming.
The reaction system is given in Table 1.1.
Table 1.1 Summary of reactions.
Reaction Number
Reaction Stoichiometry
1
xylose xylitol
2
xylitol xylulose
3
2 xylulose 3 acetaldehyde
4
acetaldehyde ethanol
5
acetaldehyde acetate
6
2 xylulose 3 glycerol
a From [1], Table 2.1, p 39.
The corresponding ODE system is [[1], p 39]
The concentrations in eqs. (1.1), denoted as [ ], are expressed in total (intraplus extracellular) moles per unit cell dry weight.
to are the kinetic rates for the six reactions listed in Table 1.1. The multiplying constants are stoichiometric coefficients. For example, reaction 3 (with rate ) in Table 1.1 produces 3 mol of acetaldehyde for every 2 mol of xylulose consumed. Therefore, eq. (1.1c) for [xylulose]/ has multiplying and eq. (1.1d) for [acetaldehyde]/ has multiplying .
The reaction rates, to , are expressed through mass action kinetics.
Note in particular the product terms for the reverse reactions in eqs. (1.2b) and (1.2c), [xylulose][ethanol] and [acetaldehyde][ethanol], which are nonlinear and therefore make the associated ODEs nonlinear (with right-hand side (RHS) terms in eqs. (1.1) that include and ). This nonlinearity precludes the usual procedures for the analytical solution of ODEs based on the linear algebra, that is, a numerical procedure is required for the solution of eqs. (1.1).
to , , , in eqs. (1.2) are kinetic constants (adjustable parameters) that are selected so that the model output matches experimental data in some manner, for example, a least squares sense. Two sets of numerical values are listed in Table 1.2
Table 1.2 Kinetic constants for two yeast strains.
BP000 refers to a wild-type yeast strain, while BP10001 refers to an engineered yeast strain.
To complete the specification of the ODE system, each of eqs. (1.1) requires an initial condition (IC) (and only one IC because these equations are first order in ).
Table 1.3 Initial conditions (ICs) for eqs. (1.1).
Equation
IC (t=0)
(1.1a)
[xylose] = 0.10724
(1.1b)
[xylitol] = 0
(1.1c)
[xylulose] = 0
(1.1d)
[acetaldehyde] = 0
(1.1e)
[ethanol] = 0
(1.1f)
[acetate] = 0
(1.1g)
[glycerol] = 0
The ODE system is now completely defined and we can proceed to programming the numerical solution.
1.3 In-Line ODE Routine
An ODE routine for eqs. (1.1) is listed in the following.
# # Library of R ODE solvers library(“deSolve”) # # Parameter values for BP10001 k1=8.87e-03; k2=13.18; k3=0.129; k4=0.497; k5=0.027; k6=0.545e-3; km2=87.7; km3=99.9; # # Initial condition yini=c(y1=0.10724,y2=0,y3=0,y4=0,y5=0,y6=0,y7=0) yini ncall=0; # # t interval nout=51 times=seq(from=0,to=2000,by=40) # # ODE programming bioreactor_1=function(t,y,parms) { with(as.list(y), { # # Assign state variables: xylose =y1; xylitol =y2; xylulose =y3; acetaldehyde=y4; ethanol =y5; acetate =y6; glycerol =y7; # # Fluxes J1=k1*xylose; J2=k2*xylitol-km2*xylulose*ethanol; J3=k3*xylulose-km3*acetaldehyde*ethanol; J4=k4*acetaldehyde; J5=k5*acetaldehyde; J6=k6*xylulose; # # Time derivatives f1=-J1; f2=J1-J2; f3=J2-2*J3-2*J6; f4=3*J3-J4-J5; f5=J4; f6=J5; f7=3*J6; # # Calls to bioreactor_1 ncall <<- ncall+1 # # Return derivative vector list(c(f1,f2,f3,f4,f5,f6,f7)) }) } # # ODE integration out=ode(y=yini,times=times,func=bioreactor_1,parms=NULL) # # ODE numerical solution for(it in 1:nout){ if(it==1){ cat(sprintf( “\n t y1 y2 y3 y4 y5 y6 y7”))} cat(sprintf(“\n %8.0f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f %8.4f”, out[it,1],out[it,2],out[it,3],out[it,4], out[it,5],out[it,6],out[it,7],out[it,8])) } # # Calls to bioreactor_1 cat(sprintf(“\n ncall = %5d\n\n”,ncall)) # # Set of 7 plots plot(out)
Listing 1.1: ODE routine.
We can note the following details about Listing 1.1.
The R library of ODE numerical integrators,
deSolve
, is specified. The contents of this library will be discussed subsequently through examples.
# # Library of R ODE solvers library(“deSolve”)
The parameters from
Table 1.2
for the engineered yeast strain
BP10001
are defined numerically.
# # Parameter values for BP10001 k1=8.87e-03; k2=13.18; k3=0.129; k4=0.497; k5=0.027; k6=0.545e-3; km2=87.7; km3=99.9;
The ICs of
Table 1.3
are defined numerically through the use of the R vector utility
c
(which defines a vector, in this case
yini
). This statement illustrates a feature of R that requires careful attention, that is, there are reserved names such as
c
that should not be used in other ways such as the definition of a variable with the name
c
.
# # Initial condition yini=c(y1=0.10724,y2=0,y3=0,y4=0,y5=0,y6=0,y7=0) yini ncall=0;
Also, the naming of the variables is open for choice (except for reserved names). Here, we select something easy to program, that is,
y1
to
y7
but programming in terms of problem-oriented variables is illustrated subsequently. Also, the elements in the IC vector
yini
are displayed by listing the name of the vector on a separate line. This is an obvious but important step to ensure that the ICs are correct as a starting point for the solution. Finally, the number of calls to the ODE function,
bioreactor_1
, is initialized.
The values of (in eqs. (1.1)) at which the solution is to be displayed are defined as the vector
times
. In this case, the R function
seq
is used to define the sequence of 51 values .
# # t interval nout=51 times=seq(from=0,to=2000,by=40)
To give good resolution (smoothness) of the plots of the solutions, 51 was selected (discussed subsequently).
Eqs. (1.1) are programmed in a function
bioreactor_1
.
# # ODE programming bioreactor_1=function(t,y,parms) { with(as.list(y), {
We can note the following details about function
bioreactor_1
.
The function is defined with three input arguments,
t,y,parms
. Also, a left brace,
{
, is used to start the function that is matched with a right brace,
}
, at the end of the function.
The input argument
y
is a list (rather than a numerical vector) specified with
with(as.list(y)
, (this statement is optional and is not used in subsequent ODE routines). The second
{
starts the
with
statement.
The seven dependent variables,
y1
to
y7
, are placed in problem-oriented variables,
xylose
to
glycerol
, to facilitate the programming of eqs. (1.1).
# # Assign state variables: xylose =y1; xylitol =y2; xylulose =y3; acetaldehyde=y4; ethanol =y5; acetate =y6; glycerol =y7;
The fluxes of eqs. (1.2) are programmed.
# # Compute fluxes J1=k1*xylose; J2=k2*xylitol-km2*xylulose*ethanol; J3=k3*xylulose-km3*acetaldehyde*ethanol; J4=k4*acetaldehyde; J5=k5*acetaldehyde; J6=k6*xylulose;
The ODEs of eqs. (1.1) are programmed, with the left-hand side (LHS) derivatives placed in the variables
f1
to
f7
. For example, [xylose]/
f1
.
# # Time derivatives f1=-J1; f2=J1-J2; f3=J2-2*J3-2*J6; f4=3*J3-J4-J5; f5=J4; f6=J5; f7=3*J6;
The number of calls to
bioreactor_1
is incremented and returned to the calling program with
<<-
.
# # Calls to bioreactor_1 ncall <<- ncall+1
This use of
<<-
illustrates a basic property of R, that is, numerical values set in a subordinate routine are not shared with higher level routines without explicit programming such as
<<-
.
The vector of derivatives is returned from bioreactor_1 as a list.
# # Return derivative vector list(c(f1,f2,f3,f4,f5,f6,f7)) }) }
Note the use of the R vector utility c. The }) ends the with statement and the second } concludes the function bioreactor_1. In other words, the derivative vector is returned from bioreactor_1 as a list. This is a requirement of the ODE integrators in the library deSolve. This completes the programming of bioreactor_1. We should note that this function is part of the program of Listing 1.1. That is, this function is in-line and is defined (programmed) before it is called (used). An alternative would be to formulate bioreactor_1 as a separate function; this is done in the next example.
Eqs. (1.1) are integrated numerically by a call to the R library integrator
ode
(which is part of
deSolve
).
# # ODE integration out=ode(y=yini,times=times,func=bioreactor_1, parms=NULL)
We can note the following details about this call to
ode
.
The inputs to
ode
are (i)
yini
, the IC vector; (ii)
times
, the vector of output values of ; and (iii)
bioreactor_1
to define the RHSs of eqs. (1.1). These inputs define the ODE system of eqs. (1.1) as expected.
The fourth input argument,
parms
, can be used to provide a vector of parameters. In the present case, it is unused. However, a vector of parameters,
k1
to
km3
, was defined previously for use in
bioreactor_1
. This sharing of the parameters with
bioreactor_1
illustrates a basic property of R: Numerical values set in a higher level routine are shared with subordinate routines (e.g., functions) without any special designation for this sharing to occur.
ode
has as a default the ODE integrator
lsoda
[10]. The
a
in the name
lsoda
stands for “automatic,” meaning that
lsoda
automatically switches between a stiff option and a nonstiff option as the numerical integration of the ODE system proceeds. The significance of stiffness will be discussed in the following and in subsequent chapters. Here we mention only that this is a sophisticated feature intended to relieve the analyst of having to specify a stiff or nonstiff integrator.
lsoda
also has a selection of options that can be specified when it is called via
ode
such as error tolerances for the ODE integration. Experimentation with these options (rather than the use of the defaults) may improve the performance of
ode
. In the present case, only the defaults are used.
The numerical solution of the ODE system is returned from
ode
as a 2D array, in this case
out
. The first index of this solution array is for the output values of the independent variable (). The second index is for the numerical solution of the ODEs. For example,
out
in the present case has the dimensions
out[51,1+7]
corresponding to (i) the 51 output values (defined previously) and (ii) the seven dependent variables of eqs. (1.1) plus the one independent variable . For example,
out[1,1]
is the value and
out[51,1]
is the value .
out[1,2]
is (from eq. (1.1a) and
Table 1.3
) [xylose] and
out[51,2]
is [xylose].
out[1,8]
is (from eq. (1.1g) and
Table 1.3
) [glycerol] and
out[51,8]
is [glycerol]. An understanding of the arrangement of the output array is essential for subsequent numerical and graphical (plotted) display of the solution.
ode
receives the number of output values of the solution from the length of the vector of output values of the independent variable. For example,
times
has 51 elements () that define the first dimension of the output array as 51 (in
out[51,1+7]
).
ode
receives the number of ODEs to be integrated from the length of the IC vector. For example,
yini
has seven elements that define the second dimension of the output array as
out[51,1+7]
(with the one added to include ).
The numerical solution is displayed at the
nout
51 output values of through a
for
loop. For
it=1
(), a heading for the numerical solution is displayed.
# # ODE numerical solution for(it in 1:nout){ if(it==1){ cat(sprintf( “\n t y1 y2 y3 y4 y5 y6 y7”))} cat(sprintf(“\n %8.0f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f %8.4f”, out[it,1],out[it,2],out[it,3],out[it,4], out[it,5],out[it,6],out[it,7],out[it,8])) }
Note the use of the values in
out
. Also, the combination of the R utilities
cat
and
sprintf
provides formatting that is used in other languages (e.g., C, C++, Matlab).
The number of calls to
bioreactor_1
is displayed at the end of the solution to give an indication of the computational effort required to compute the solution.
# # Calls to bioreactor_1 cat(sprintf(“\n ncall = %5d\n\n”,ncall))
Finally, the solutions of eqs. (1.1) are plotted with the R utility
plot
.
# # Set of 7 plots plot(out)
A complete plot is produced with just this abbreviated use of
out
.
plot
has a variety of options to format the graphical output that will be considered in subsequent applications.
1.4 Numerical and Graphical Outputs
Abbreviated numerical output from Listing 1.1 is given in Table 1.4.
Table 1.4 Abbreviated numerical output from Listing 1.1.
t y1 y2 y3 y4 y5 y6 y7 0 0.1072 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 40 0.0752 0.0020 0.0153 0.0009 0.0195 0.0011 0.0006 80 0.0527 0.0053 0.0221 0.0008 0.0361 0.0020 0.0018 120 0.0370 0.0081 0.0245 0.0006 0.0497 0.0027 0.0034 160 0.0259 0.0101 0.0248 0.0005 0.0608 0.0033 0.0050 200 0.0182 0.0112 0.0239 0.0004 0.0702 0.0038 0.0066 . . . . . . Output for t = 240 to 1760 removed . . . . . . 1800 0.0000 0.0004 0.0004 0.0000 0.1303 0.0071 0.0222 1840 0.0000 0.0003 0.0004 0.0000 0.1304 0.0071 0.0223 1880 0.0000 0.0003 0.0004 0.0000 0.1305 0.0071 0.0223 1920 0.0000 0.0003 0.0003 0.0000 0.1306 0.0071 0.0223 1960 0.0000 0.0003 0.0003 0.0000 0.1306 0.0071 0.0223 2000 0.0000 0.0002 0.0003 0.0000 0.1307 0.0071 0.0223 ncall = 427
We can note the following details about this output.
The ICs (at ) correspond to the values in
Table 1.3
. While this may seem to be an obvious fact, it is a worthwhile check to ensure that the solution has the correct starting values.
The solutions approach steady-state conditions as . Note in particular that
y1
(for xylose from eq. (1.1a)) approaches zero as the the reactant that drives the system is nearly consumed. Also,
y5
(for ethanol from eq. (1.1e)) approaches 0.1307 indicating a significant production of ethanol, the product of primary interest (e.g., possibly to be used as a fuel).
y7
(for glycerol from eq. (1.1g)) approaches 0.0223 and might represent a contaminant that would have to be subsequently reduced by a separation process; this is rather typical of reaction systems, that is, they usually produce undesirable by-products.
The computational effort is quite modest,
ncall = 427
(the reason for calling this “modest” is explained subsequently).
The graphical output is given in Fig. 1.1. We can note the following about Fig. 1.1.
The plotting utility
plot
provides automatic scaling of each of the seven dependent variables. Also, the default of
plot
is the solid lines connecting the values in
Table 1.4
; alternative options provide discrete points, or points connected by lines.
The initial () values reflect the ICs of
Table 1.3
and the final values () reflect the values of
Table 1.4
.
The solutions have their largest derivatives at the beginning which is typical of ODE systems (the LHSs of eqs. (1.2) is largest initially).
The plots are smooth with 51 points.
Figure 1.1 Solutions to eqs. (1.1).
A fundamental question remains concerning the accuracy of the solution in Table 1.4. As an exact (i.e., analytical, mathematical, closed form) solution is not available for eqs. (1.1) (primarily because they are nonlinear as discussed previously), we cannotdirectly determine the accuracy of the numerical solution by comparison with an exact solution (and if such a solution was available, there would really be no need to compute a numerical solution).
We therefore must use a method of accuracy evaluation that is built on the numerical approach. For example, we could change the specified error tolerances for lsoda (via the call to ode) and compare the solutions as the error tolerances are changed. Or we could use other ODE integrators (other than lsoda) and compare the solutions from different integrators (this approach is discussed in a subsequent example). In any case, some form of error analysis is an essential part of any numerical procedure to give reasonable confidence that the numerical solution has acceptable accuracy.
Finally, with the operational code of Listing 1.1, we can now perform studies (experiments) that will contribute to an understanding of the problem system (which is usually the ultimate objective in developing a mathematical model) on the computer. For example, the effect of changing the model parameters, termed the parameter sensitivity, can be carried out by observing the changes in the solutions as parameters are varied. As an example, the BP000 parameters of Table 1.2 can be used in place of the BP10001 parameters (in Listing 1.1) to investigate the effect of using an engineered yeast strain (BP10001) in place of a wild-type yeast strain (BP000). Ideally, an increase in ethanol production would be observed (the final value of y5 in Table 1.4 would increase), indicating that the engineered yeast strain can improve the efficiency of ethanol production.
This type of parameter sensitivity analysis presupposes available values of the model parameters that reflect the performance of the problem system, and these parameters might have to be measured experimentally, for example, by comparing the model solution with laboratory data, and/or estimated using available theory. A good example of the comparison of the ethanol model solution with experimental data is given in [1] for BP000 (Fig. 2.7) and BP10001 (Fig. 2.8).
1.5 Separate ODE Routine
Variations of the coding in Listing 1.1 will now be considered. The intent is to produce a more flexible modular format and to enhance the graphical output. The main program now is in Listing 1.2 (in place of Listing 1.1)
# # Library of R ODE solvers library(“deSolve”) # # ODE routine setwd(“c:/R/bme_ode/chap1”) source(“bioreactor_2.R”) # # Parameter values for BP10001 k1=8.87e-03; k2=13.18; k3=0.129; k4=0.497; k5=0.027; k6=0.545e-3; km2=87.7; km3=99.9; # # Initial condition yini=c(0.10724,0,0,0,0,0,0) ncall=0; # # t interval nout=51 times=seq(from=0,to=2000,by=40) # # ODE integration out=ode(y=yini,times=times,func=bioreactor_2,parms=NULL) # # ODE numerical solution for(it in 1:nout){ if(it==1){ cat(sprintf( “\n t y1 y2 y3 y4 y5 y6 y7”))} cat(sprintf(“\n %8.0f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f %8.4f”, out[it,1],out[it,2],out[it,3],out[it,4], out[it,5],out[it,6],out[it,7],out[it,8])) } # # Calls to bioreactor_2 cat(sprintf(“\n ncall = %5d\n\n”,ncall)) # # Single plot par(mfrow=c(1,1)) # # y1 plot(out[,1],out[,2],type="l",xlab="t",ylab="y1(t),..., y7(t)", xlim=c(0,2000),ylim=c(0,0.14),lty=1, main="y1(t),..., y7(t) vs t", lwd=2) # # y2 lines(out[,1],out[,3],type="l",lty=2,lwd=2) # # y3 lines(out[,1],out[,4],type="l",lty=3,lwd=2) # # y4 lines(out[,1],out[,5],type="l",lty=4,lwd=2) # # y5 lines(out[,1],out[,6],type="l",lty=5,lwd=2) # # y6 lines(out[,1],out[,7],type="l",lty=6,lwd=2) # # y7 lines(out[,1],out[,8],type="l",lty=7,lwd=2)
Listing 1.2: Main program with separate ODE routine.
We can note the following details about Listing 1.2.
library(“deSolve”)
is used again (as in Listing 1.1) in order to access the ODE integrator
ode
. In addition, the separate ODE routine
bioreactor_2
is accessed through the
setwd
and
source
R utilities.
# # Library of R ODE solvers library(“deSolve”) # # ODE routine setwd(“c:/R/bme_ode/chap1”) source(“bioreactor_2.R”)
To explain the use of
setwd
and
source
:
setwd, set woking directory, is used to go to a directory (folder) where the R routines are located. Note in particular the use of the forward slash / rather than the usual backslash .
source
identifies a particular file within the directory identified by the
setwd
; in this case, the ODE routine
bioreactor_2
is called by
ode
.
These two statements could be combined as
source(“c:/R/bme_ode/chap1/bioreactor_2.R”)
If the R application uses a series of files from the same directory, using the
setwd
is usually simpler; a series of
source
statements can then be used access the required files.
The sections of Listing 1.2 for setting the parameters, IC and interval, are the same as in Listing 1.1 and are therefore not discussed here. The call to
ode
uses the ODE routine
bioreactor_2
(Listing 1.3) rather than
bioreactor_1
(Listing 1.1).
# # ODE integration out=ode(y=yini,times=times,func=bioreactor_2, parms=NULL)
Again, the ODE solution is returned in 2D array
out
for subsequent display.
bioreactor_2
is in a separate routine rather than placed in-line as in Listing 1.1, which makes the coding more modular and easier to follow.
The display of the numerical solution is the same as in Listing 1.1 so this code is not discussed here.
The number of calls to
bioreactor_3
(returned from
bioreactor_2
at the end of the solution, i.e., at ) is displayed.
# # Calls to bioreactor_2 cat(sprintf(“\n ncall = %5d\n\n”,ncall))
The graphical output is extended to produce a single plot with the seven ODE solution curves.
# # Single plot par(mfrow=c(1,1)) # # y1 plot(out[,1],out[,2],type="l",xlab="t",ylab="y1(t), ...,y7(t)", xlim=c(0,2000),ylim=c(0,0.14),lty=1, main="y1(t), ...,y7(t) vs t", lwd=2) # # y2 lines(out[,1],out[,3],type="l",lty=2,lwd=2) # # y3 lines(out[,1],out[,4],type="l",lty=3,lwd=2) # # y4 lines(out[,1],out[,5],type="l",lty=4,lwd=2) # # y5 lines(out[,1],out[,6],type="l",lty=5,lwd=2) # # y6 lines(out[,1],out[,7],type="l",lty=6,lwd=2) # # y7 lines(out[,1],out[,8],type="l",lty=7,lwd=2)
To explain this coding,
A array of plots is specified, that is, a single plot;
# # Single plot par(mfrow=c(1,1))
plot
is used with a series of parameters for .
# # y1 plot(out[,1],out[,2],type="l",xlab="t",ylab="y1 (t),...,y7(t)", xlim=c(0,2000),ylim=c(0,0.14),lty=1, main="y1 (t),...,y7(t) vs t", lwd=2)
These parameters are:
out[,1],out[,2]
plotted to give a solution curve for eq. (1.1a) of versus ;
type=“l”
designates a line type of the solution curve (rather than a point type);
xlab=“t”
specifies the label
t
on the abcissa (horizontal) axis;
ylab=“y1(t),…,y7(t)”
specifies the label on the ordinate (vertical) axis;
xlim=c(0,2000)
scales the horizontal axis for ;
ylim=c(0,0.14)
scales the vertical axis to include the range of values from to ;
lty=1
sets the type of line for the first solution as reflected in
Fig. 1.2
;
main=“y1(t),...,y7(t) vs t”
specifies a main label or title for the plot as reflected in
Fig. 1.2
;
lwd=2
sets the line width for the first solution as reflected in
Fig. 1.2
.
is included as a second solution with the R utility
lines
by plotting
out[,1],out[,3]
. The parameters are the same as for the previous call to
plot
except
lty=2
, which specifies a second type of line as reflected in
Fig. 1.2
.
# # y2 lines(out[,1],out[,3],type="l",lty=2,lwd=2)
to are plotted in the same way with
lines
. For example, is plotted as
# # y7 lines(out[,1],out[,8],type="l",lty=7,lwd=2)
with a line type specified as
lty=7
, which specifies a seventh type of line as reflected in
Fig. 1.2
.
Figure 1.2 Solutions to eqs. (1.1) using a separate ODE routine.
bioreactor_2 in Listing 1.3 is a separate routine called by ode.
bioreactor_2=function(t,y,parms) { # # Assign state variables: xylose =y[1]; xylitol =y[2]; xylulose =y[3]; acetaldehyde=y[4]; ethanol =y[5]; acetate =y[6]; glycerol =y[7]; # # Compute fluxes J1=k1*xylose; J2=k2*xylitol-km2*xylulose*ethanol; J3=k3*xylulose-km3*acetaldehyde*ethanol; J4=k4*acetaldehyde; J5=k5*acetaldehyde; J6=k6*xylulose; # # Time derivatives f1=-J1; f2=J1-J2; f3=J2-2*J3-2*J6; f4=3*J3-J4-J5; f5=J4; f6=J5; f7=3*J6; # # Calls to bioreactor_2 ncall <<- ncall+1 # # Return derivative vector return(list(c(f1,f2,f3,f4,f5,f6,f7))) }
Listing 1.3: ODE routine bioreactor_2.
bioreactor_2 is the same as bioreactor_1 of Listing 1.1 except for the following details.
The function is defined as in Listing 1.1, but the statement specifying
y
as a list (
with(as.list(y))
is not used.
bioreactor_2=function(t,y,parms) {
The dependent variables constitute a vector (
y[1],...,y[7]
) rather than a list of scalars (
y1,...,y7
) as in Listing 1.1. In other words, the input argument of
bioreactor_2
,
y
, is a vector and not a list.
At the end, the calls to
bioreactor_3
is incremented and returned to the main program of Listing 1.2 with
<<-
.
# # Calls to bioreactor_2 ncall <<- ncall+1
Finally, the derivatives
f1
to
f7
are placed in a vector with
c
that is returned from
bioreactor_2
(to
lsoda
in
ode
) as a list (as required by the ODE integrators in
deSolve
).
# # Return derivative vector return(list(c(f1,f2,f3,f4,f5,f6,f7))) }
The right-hand bracket
}
concludes
bioreactor_2
.
The numerical output from Listing 1.2 is identical to that of Listing 1.1 as expected because the only difference in the coding is the use of the separate ODE routine bioreactor_2. The graphical output is given in Fig. 1.2. Note the composite plot of Fig 1.2 rather than the separate plots of Fig 1.1.
This concludes the example with a separate ODE routine. The intent is primarily to demonstrate the use of subordinate functions to modularize the code (rather than having it all in one routine such as in Listing 1.1). This modularization becomes increasingly useful as the complexity of the application increases because it can be used to organize the code into small, more easily manageable sections.
1.6 Alternative Forms of ODE Coding
We now consider variations of the coding in the preceding Listings 1.1–1.3. The intention is primarily to introduce alternatives that can be useful, particularly as the number of ODEs increases (beyond the 7 of eqs. (1.1)).
In the following example, the main program is the same as in Listing 1.2 except for the use of bioreactor_3 in place of bioreactor_2.
# # ODE integration out=ode(y=yini,times=times,func=bioreactor_3,parms=NULL)
bioreactor_3 is similar to bioreactor_2 in Listing 1.3, except that the fluxes and the derivatives are programmed as vectors (Listing 1.4).
bioreactor_3=function(t,y,parms) { # # Assign state variables: xylose =y[1]; xylitol =y[2]; xylulose =y[3]; acetaldehyde=y[4]; ethanol =y[5]; acetate =y[6]; glycerol =y[7]; # # Compute fluxes J=rep(0,n) J[1]=k1*xylose; J[2]=k2*xylitol-km2*xylulose*ethanol; J[3]=k3*xylulose-km3*acetaldehyde*ethanol; J[4]=k4*acetaldehyde; J[5]=k5*acetaldehyde; J[6]=k6*xylulose; # # Time derivatives f=rep(0,n) f[1]=-J[1]; f[2]=J[1]-J[2]; f[3]=J[2]-2*J[3]-2*J[6]; f[4]=3*J[3]-J[4]-J[5]; f[5]=J[4]; f[6]=J[5]; f[7]=3*J[6]; # # Calls to bioreactor_3 ncall <<- ncall+1 # # Return derivative vector return(list(c(f))) }
Listing 1.4: ODE routine with vectors to facilitate the ODE programming.
We can note the following details about bioreactor_3.
The sizes of the vectors for the fluxes and derivatives are declared using the R utility
rep
before the vectors are used.
