Chapter 1
Introduction 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.

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).

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.

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...