King M.R., Mody N.A. Numerical and Statistical Methods for Bioengineering: Applications in MATLAB

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

of drug absorbed by the body, irrespective of the rate of absorption. This is useful

when trying to determine whether two formulations of the same dose (for example a

capsule and a tablet) release the same dose of drug to the body. Another use is in the

therapeutic monitoring of toxic drugs. For example, gentamicin is an antibiotic which

displays nephro- and ototoxicities; measurement of gentamicin concentrations in a

patient’s plasma, and calculation of the AUC is used to guide the dosage of this drug.

(a) Suppose you want to design an experiment to calculate the AUC in a patient

over a four-hour period, using a four-point Gaussian quadrature. At what times

(starting the experiment at t = 0) should you make the plasma measurements?

Show your work.

(b) Now redesign the experiment to obtain plasma measurements appropriate for

using an eight-panel Simpson’s rule quadrature.

Weights and nodes for the Gaussian four-point rule are given in Table P6.5.

6.9. Reformulate Equation (6.23) to show how the trapezoidal formula for step size h/2,

Ih=2ðÞ, is related to the trapezoidal formula for step size h. Rewrite Program 6.5 to

use your new formula to generate the ﬁrst column of approximations for Romberg

integration. You will not need to call the user-deﬁned MATLAB function

trapezoidal_rule since your program will itself perform the trapezoidal inte-

gration using the improvised formula. Verify the correctness of your program by

using it to calculate the time taken for the IV drip bag to empty by 90% (see Box 6.3).

6.10. Rewrite the function Romberg_integration in Program 6.5 so that levels of

integration are progressively increased until the speciﬁed tolerance is met according

to the stopping criterion

j;k

 R

j;k1

j;k

 Tol:

The tolerance level should be an input parameter for the function. Thus, the

program will determine the maximum level required to obtain a solution with the

desired accuracy. You will need to modify the program so that an additional row of

integral approximations corresponding to the halved step size (including the column

1 trapezoidal rule derived approximation) is dynamically generated when the result

does not meet the tolerance criterion. You may use the trapezoidal_rule

function listed in Program 6.1 to perform the numerical integration step, or the

new formula you derived in Problem 6.9 to generate the ﬁrst column of approxima-

tions. Repeat the IV drip calculations discussed in Box 6.3 for tolerances Tol = 10

−7

−8

, and 10

−9

6.11. Rewrite Program 6.6 to perform an (n + 1)-point Gauss–Legendre quadrature,

where n = 1, 2, 3, or 4. Rework Example 6.3 for n = 2, 3, and 4 and N (number

of segments) = 2. You may wish to use the switch-case program control state-

ment to assign the appropriate nodes and weights to the node and weight variable

Table P6.5. Weight s and nodes for the

Gaussian four-point rule

±0.8611363 0.3478548

±0.3399810 0.6521452

407

6.7 Problems

based on user input. See Appendix A or MATLAB help for the appropriate syntax

for switch-case. This control statement should be used in place of the if-else-

end statement when several different outcomes are possible, depending on the value

of the control variable.

References

Abramowitz, M. and Stegun, I. A. (1965) Handbook of Mathematical Functions (Mineola,

NY: Dover Publications, Inc.).

Anderson, J. L. and Quinn, J. A. (1974) Restricted Transport in Small Pores. A Model for

Steric Exclusion and Hindered Particle Motion. Biophys. J., 14, 130–50.

Bailey, P. L., Lu, J. K., Pace, N. L. et al. (2000) Effects of Intrathecal Morphine on the

Ventilatory Response to Hypoxia. N. Engl. J. Med., 343, 1228–34.

Brenner, H. and Gaydos, L. J. (1977) Constrained Brownian-Movement of Spherical-

Particles in Cylindrical Pores of Comparable Radius – Models of Diffusive and

Convective Transport of Solute Molecules in Membranes and Porous-Media. J.

Colloid Interface Sci., 58, 312–56.

Burden, R. L. and Faires, J. D. (2005) Numerical Analysis (Belmont, CA: Thomson Brooks/

Cole).

Deen, W. M. (1987) Hindered Transport of Large Molecules in Liquid-Filled Pores. Aiche

J., 33, 1409–25.

Grasselli, M. and Pelinovsky, D. (2008) Numerical Mathematics (Sudbury, MA: Jones and

Bartlett Publishers).

Patel, V. A. (1994) Numerical Analysis (Fort Worth, TX: Saunders College Publishing).

Ralston, A. and Rabinowitz, P. (1978) A First Course in Numerical Analysis (New York:

McGraw-Hill, Inc.).

Recktenwald, G. W. (2000) Numerical Methods with Matlab Implementation and Application

(Upper Saddle River, NJ: Prentice Hall).

Richter, W. and Richter, M. (2003) The Shape of the FMRI Bold Response in Children and

Adults Changes Systematically with Age. Neuroimage, 20, 1122–31.

Smith, F. G. and Deen, W. M. (1983) Electrostatic Effects on the Partitioning of Spherical

Colloids between Dilute Bulk Solution and Cylindrical Pores. J. Colloid Interface Sci.,

91, 571–90.

408

Numerical quadrature

7 Numerical integration of ordinary

differential equations

7.1 Introduction

Modeling of dynamic processes, i.e. the behavior of physical, chemical, or bio-

logical systems that have not reached equilibrium, is done using mathematical

formulations that contain one or more differential terms such as dy=dx or dC=dt.

Such equations are called differential equations. The derivative in the equation

indicates that a quantity is constantly evolving or changing with respect to time t or

space (x, y,orz). If the quantity of interest varies with respect to a single variable

only, then the equation has only one independent variable and is called an ordinary

differential equation (ODE). An ODE can have more than one dependent variable,

but all dep endent variables in the equation vary as a function of the same inde-

pendent variable. In other words, all derivatives in the equation are taken with

respect to one independent variable. An ODE can be contrasted with a partial

differential equation (PDE). A PDE contains two or more partial derivatives since

the dependent variable is a function of more than one independent variable. An

example of an ODE and a PDE are given below. We will not concern ourselves with

PDEs in this chapter.

ODE:

¼ 3 þ 4y;

PDE:

∂T

∂t

¼ 10

∂

∂x

The order of an ordinary differential equation is equal to the order of the highest

derivative contained in the equation. A differential equation of ﬁrst order con-

tains only the ﬁrst derivative of the dependent variable y, i.e. dy=dt.Annth-order

ODE contains the nth derivative, i.e. d

y=dt

. All ODEs for which n ≥ 2 are

collectively termed higher-order ordinary differential equations. An analytical sol-

ution of an ODE describes the dependent variable y as a function of the

independent variable t such that the functional form y ¼ fðtÞ satisﬁes the differ-

ential equation. A unique solution of the ODE is implied only if n constraints or

boundary conditions are deﬁned along with the equation, where n is equal to the

order of the ODE. Why is this so? You know that the analytical solution of a

ﬁrst-order differential equation is obtained by performing a single integration

step, which produces one integration constant. Solving an nth-order differential

equation is synonymous with performing n sequential integrations or producing n

constants. To determine the values of the n constants, we need to specify n

conditions.

ODEs for which all n co nstraints apply at the beginning of the integration interval

are solved as initial-value problems. Most initial-value ODE problems have time-

dependent variables. An ODE that contains a time derivative represents a process

that has not reached steady state or equilibrium.

Many biological and industrial

processes are time-dependent. We use the t symbol to represent the independent

variable of an initial-value problem (IVP). A ﬁrst-order ODE is always solved as

an IVP since only one boundary condition is required, which provides the initial

condition, y(0) = y

. An initial-value problem consisting of a higher-order ODE,

such as

þ xy

;

requires that, in addition to the initial condition y 0ðÞ¼y

, initial values of the

derivatives at t = 0, such as y

0ðÞ¼a; y

0ðÞ¼b; must also be speciﬁed to obtain a

solution. Sometimes boundary conditions are speciﬁed at two points of the problem,

i.e. at the beginning and ending time points or spatial positions. Such problems are

called two-point boundary-value problems. A different solution method must be used

for such problems.

When the differential equation exhibits linear dependence on the derivative terms

and on the dependent variable, it is called a linear ODE. All derivatives and terms

involving the dependent variable in a linear equation have a power of 1. An example

of a linear ODE is

þ 3

þ t

sin t þ ty ¼ 0:

Note that there is no requirement for a linear ODE to be linearly dependent on the

independent variable t. If some of the dependent variable terms are multiplied with

each other, or are raised to multiple powers, or if the equation has transcedental

functions of y, e.g. sin(y) or exp(y), the resulting ODE is nonlinear.

Analytical techniques are available to solve linear and nonlinear differential

equations. However, not all differential equations are amenable to an analytical

solution, and in practice this is more often the case. Sometimes an ODE may contain

more than one depend ent variable. In these situations, the ODE must be coupled

with other algebraic equations or ODEs before a solution can be found. Several

coupled ODEs produce a system of simultaneous ODEs. For example,

¼ 3y þ z þ 1;

¼ 4z þ 3;

form a set of two simultaneous ﬁrst-order linear ODEs. Finding an analytical

solution to coupled differential equations can be even more challenging. Instead of

A steady state process can involve a spatial gradient of a quantity, but the gradient is unchanging with

time. The net driving forces are ﬁnite and maintain the gradient. There is no accumulation of depletion

of the quantity of interest anywhere in the system. On the other hand, a system at equilibrium has no

spatial gradient of any quantity. The net driving forces are zero.

410

Numerical integration of ODEs

searching for an elusive analytical solution, it can be more efﬁcient to ﬁnd the unique

solution of a well-behaved (well-posed) ODE problem using numerical techniques.

Any numerical method of solution for an initial-value or a boundary-value problem

involves discr etizing the continuous time domain (or space) into a set of ﬁnite

discrete points (equally spaced or unequally spaced) or nodes at which we ﬁnd the

magnitude of yðtÞ. The value of y evolves in time according to the slope of the curve

prescribed by the right-hand side of the ODE, dy=dt ¼ fðt; yÞ. The solution at each

node is obtained by marching step-wise forward in time. We begin where the solution

is known, such as at t = 0, and use the known solutions at one or more nearby points

to calculate the unknown solution at the adjacent point. If a solution is desired at an

intermediate point that is not a node, it can be obtained using interpolation. The

numerical solution is only approximate because it suffers from an accumulation of

truncation errors and round-off errors.

In this chapter, we will discuss several classical and widely used numerical

methods such as Euler methods, Runge–Kutta methods, and the predictor–corrector

methods based on the Adam s–Bashforth and Adams–Moulton methods for solving

ordinary differential equations. For each of these methods, we will investigate the

truncation error. The order of the truncation error determines the convergence rate

of the error to zero and the accuracy limitations of the numerical method. Even when

a well-behaved ODE has a bounded solution, a numerical solver may not necessarily

ﬁnd it. The conditions, or values of the step size for which a numerical method is

stable, i.e. conditions under which a bounded solution can be found, are known as

numerical stability criteria. These depend upon both the nature of the differential

problem and the characteristics of the numerical solution method. We will consider

the stability issues of each method that we discuss.

Because we are focused on introducing and applying numerical methods to solve

engineering problems, we do not dwell on the rigors of the mathematics such as

theorems of uniqueness and existence and their corresponding proofs. In-depth

mathematical analyses of numerical methods are offered by textbooks such as

Numerical Analysis by Burden and Faires (2005), and Numerical Mathematics by

Grasselli and Pelinovsky (2008). In our discussion of numerical methods for ODEs,

it is assumed (unless mentioned otherwise) that the dependent variable yðtÞ is

continuous and fully differentiable within the domain of interest, the derivatives of

yðtÞ of successive orders y

tðÞ; y

tðÞ; ... are continuous within the deﬁned domain,

and that a unique bounded solution of the ODE or ODE system exists for the

boundary (initial) conditions speciﬁed.

In this chapter, we assume that the round-off error is much smaller than the

truncation error. Accordingly, we neglect the former as a component of the overall

error. When numerical precision is less than double, round-off error can sizably

contribute to the total error, decreasing the accuracy of the solution. Since

MATLAB uses double precision by default and all our calculations will be per-

formed using MATLAB software, we are justiﬁed in neglecting round-off error. In

Sections 7.2–7.5, methods to solve initial-value problems consisting of a single

equation or a set of coupled equations are discussed in detail. The numerical method

chosen to solve an ODE problem will depend on the inherent properties of the

equation and the limitations of the ODE solver. How one should approach stiff

ODEs, which are more challenging to solve, is the topic of discussion in Section 7.6.

Section 7.7 covers boundary-value problems and focuses on the shooting method, a

solution technique that converts a boundary-value problem to an initial-value

problem and uses iterative techniques to ﬁnd the solution.

411

7.1 Introduction

Box 7.1A HIV–1 dynamics in the blood stream

HIV–1 virus particles attack lymphocytes and hijack their genetic machinery to manufacture many virus

particles within the cell. Once the viral particles have been prepared and packaged, the viral genome

programs the infected host to undergo lysis. The newly assembled virions are released into the blood

stream and are capable of attacking new cells and spreading the infection. The dynamics of viral

infection and replication in plasma can be modeled by a set of differential equations (Perelson et al.,

1996). If T is the concentration of target cells, T



is the concentration of infected cells, and V is the

concentration of infectious virus particles in plasma, we can write



¼ kVT  δT



;

¼ NδT



 cV;

where k characterizes the rate at which the target cells T are infected by viral particles, δ is the rate of cell

loss by lysis, apoptosis, or removal by the immune system, c is the rate at which viral particles are

cleared from the blood stream, and N is the number of virions produced by an infected cell.

In a clinical study, an HIV–1 protease inhibitor, ritonavir, was administered to five HIV-infected

individuals (Perelson et al., 1996). Once the drug is absorbed into the body and into the target cells, the

newly produced virus particles lose their ability to infect cells. Thus, once the drug takes effect, all newly

produced virus particles are no longer infectious. The drug does not stop infectious viral particles from

entering into target cells, or prevent the production of virus by existing infected cells. The dynamics of

ritonavir-induced loss of infectious virus particles in plasma can be modeled with the following set of

first-order ODEs:



¼ kV

T  δT



; (7:1)

¼cV

; (7:2)

¼ NδT



 cV

; (7:3)

where V

is the concentration of infectious viral RNA in plasma, and V

is the concentration of non-

infectious viral particles in plasma. The experimental study yielded the following estimates of the

reaction rate parameters (Perelson et al., 1996):

δ ¼ 0:5=day;

c ¼ 3:0=day:

The initial concentrations of T, T



, V

, and V

are as follows:

t ¼ 0ðÞ¼100/μl;

t ¼ 0ðÞ¼0/μl;

Tt¼ 0ðÞ¼250 non-infected cells/μl (Haase et al., 1996);



t ¼ 0ðÞ¼10 infected cells/ μl (Haase et al., 1996).

Based on a quasi-steady state analysis (dT



=dt ¼ 0, dV=dt ¼ 0) before time t = 0, we calculate

the following:

k ¼ 2  10

4

μl/day/virions, and

N = 60 virions produced per cell.

Using these equations, we will investigate the dynamics of HIV infection following protease inhibitor

therapy.

412

Numerical integration of ODEs

Box 7.2A Enzyme deactivation in industrial processes

Enzymes are proteins that catalyze physiological reactions, typically at body temperatur e, atmospheric

pressur e, and within a very specific pH range. The process industry has made significant advances in enzyme

catalysis technology. Adoption of biological processes for industrial purposes has several key advantages over

using traditional chemical methods for manufacture: (1) inherently safe operating conditions are characteristic

of enzyme-catalyzed reactions, and (2) non-toxic non-polluting waste products are produced by enzymatic

react ions. Applications of enzymes include manufacture of food, drugs, and deter gents, carrying out

intermediate steps in synthesis of chemicals, and processing of meat, hide, and alcoholic beverages.

Enzymes used in an industrial setting begin to deactivate over time, i.e. lose their activity. Enzyme

denaturation may occur due to elevated process temperatures, mechanically abrasive operations that

generate shearing stresses on the molecules, solution pH, and/or chemical denaturants. Sometimes

trace amounts of a poison present in the broth or protein solution may bind to the active site of the

enzymes and slowly but irreversibly inactivate the entire batch of enzymes.

A simple enzyme-catalyzed reaction can be modeled by the following reaction sequence:

E þ S Ð ES !

E þ P;

E: uncomplexed or free active enzyme,

S: substrate,

ES: enzyme–substrate complex,

P: product.

The reaction rate is given by the Michaelis–Menten equation,



¼ r ¼

total

þ S

; (7:4)

where E

total

is the total concentration of active enzyme, i.e. E

total

¼ E þ ES.

The concentration of active enzyme progressively reduces due to gradual deactivation. Suppose that

the deactivation kinetics follow a first-order process:

total

¼k

It is assumed that only the free (unbound) enzyme can deactivate. We can write E ¼ E

total

 ES, and we

have from Michaelis–Menten theory

ES ¼

total

þ S

(An important assumption made by Michaelis and Menten was that dðESÞ=dt ¼ 0. This is a valid

assumption only when the substrate concentration is much greater than the total enzyme concentration.)

The above kinetic equation for dE

total

=dt becomes

total

¼

total

1 þ S=K

: (7:5)

Equations (7.4) and (7.5) are two coupled first-order ordinary differential equations. Pr oduct formation and

enzyme deactivation are functions of the substrate concentration and the total active enzyme concentration.

An enzyme catalyzes proteolysis (hydrolysis) of a substrate for which the rate parameters are

¼ 21 M/s/e.u. (e.u. ≡ enzyme units);

¼ 0:004 M;

¼ 0:03/s.

The starting substrate and enzyme concentrations are S

¼ 0:05 M and E

¼ 1  10

6

units. How

long will it take for the total active enzyme concentration to fall to half its initial value? Plot the change in

substrate and active enzyme concentrations with time.

413

7.1 Introduction

Box 7.3A Dynamics of an epidemic outbreak

When considering the outbreak of a communicable disease (e.g. swine flu), members of the affected

population can be categorized into three groups: diseased individuals D, susceptible individuals S, and

immune individuals I. In this particular epidemic model, individuals with immunity contracted the

disease and recover from it, and cannot become reinfected or transmit the disease to another person.

Depending on the mechanism and rate at which the disease spreads, the disease either washes out, or

eventually spreads throughout the entire population, i.e. the disease reaches epidemic proportions. Let k

be a rate quantity that characterizes the ease with which the infection spreads from an infected person to

an uninfected person. Let the birth rate per individual within the population be β and the death rate be μ.

The change in the number of susceptible individuals in the population can be modeled by the

deterministic equation

¼kSD þ βðS þ IÞμS: (7:6)

Newborns are always born into the susceptible population. When a person recovers from the disease, he

or she is permanently immune. A person remains sick for a time period of γ, which is the average

infectious period.

The dynamics of infection are represented by the following equations:

¼ kSD 

D; (7:7)

D  μI: (7:8)

It is assumed that the recovery period is short enough that one does not generally die of natural causes

while infected. We will study the behavior of these equations for two sets of values of the defining

parameters of the epidemic model.

Box 7.4A Microbial population dynamics

In ecological systems, multiple species interact with each other in a variety of ways. When animals compete

for the same food, this interaction is called competition. A carnivorous animal may prey on a herbivore in a

predator–prey type of interaction. Several distinct forms of interaction are observed among microbial

mixed-culture populations. Two organisms may exchange essential nutrients to benefit one another mutually.

When each species of the mixed population fares better in the presence of the other species as opposed to

under isolated conditions, the population interaction dynamics is called mutualism. Sometimes, in a two-

population system, only one of the two species may derive an advantage (commensalism). For example, a

yeast may metabolize glucose to pr oduce an end product that may serve as bacterial nourishment. The yeast

however may not derive any benefit from co-existing with the bacterial species. In another form of interaction,

the presence of one species may harm or inhibit the growth of another, e.g. the metabolite of one species

may be toxic to the other. In this case, growth of one species in the mixed culture is adversely affected, while

the growth of the second species is unaffected (amensalism).

Microbes play a variety of indispensable roles both in the natural environment and in industrial

processes. Bacteria and fungi decompose organic matter such as cellulose and proteins into simpler

carbon compounds that are released to the atmosphere in inorganic form as CO

and bicarbonates. They

are responsible for recirculating N, S, and O through the atmosphere, for instance by fixing atmospheric

nitrogen, and by breaking down proteins and nucleic acids into simpler nitrogen compounds such as

, nitrites, and nitrates. Mixed cultures of microbes are used in cheese production, in fermentation

broths to manufacture alcoholic beverages such as whiskey, rum, and beer, and to manufacture

vitamins. Domestic (sewage) wastewater and industrial aqueous waste must be treated before being

released into natural water bodies. Processes such as activated sludge, trickling biological filtration, and

414

Numerical integration of ODEs

anaerobic digestion of sludge make use of a complex mix of sludge, sewage bacteria, and bacteria-

consuming protozoa, to digest the organic and inorganic waste.

A well-known mathematical model that simulates the oscillatory behavior of predator and prey

populations (e.g. plant –herbivore, insect pest–control agent, and herbivore–carnivore) in a two-species

ecological system is the Lotka–Volterra model.LetN

be the number of prey in an isolated habitat,

and let N

be the number of predators localized to the same region. We make the assumption that the

prey is the only source of food for the predator. According to this model, the birth rate of prey is

proportional to the size of its population, while the consumption rate of prey is proportional to the

number of predators multiplied by the number of prey. The Lotka–Volterra model for predator–prey

dynamics is presented below:

¼ bN

 γN

; (7:9)

¼ εγN

 dN

: (7:10)

The parameter γ describes the ease or difficulty with which a prey can be killed and devoured by the

predator, and together the term γN

equals the frequency of killing events. The survival and growth of

the predator population depends on the quantity of food available (N

) and the size of the predator

population (N

); ε is a measure of how quickly the predator population reproduces for each successfully

hunted prey; d is the death rate for the predator population. A drawback of this model is that it assumes

an unlimited supply of food for the prey. Thus, in the absence of predators, the prey population grows

exponentially with no bounds.

A modified version of the Lotka–Volterra model uses the Monod equation (see Box 1.6 for a

discussion on Monod growth kinetics) to model the growth of bacterial prey as a nonlinear function

of substrate concentration, and the growth of a microbial predator as a function of prey concentration

(Bailey and Ollis, 1986). This model was developed by Tsuchiya et al.(1972) for mixed population

microbial dynamics in a chemostat (a continuous well-stirred flow reactor containing a uniform

suspension of microbes). A mass balance performed on the substrate over the chemostat apparatus

(see Figure P7.2 in Problem 7.5), assuming constant density and equal flowrates for the feed stream

and product stream, yields the following time-dependent differential equation:

¼ FS

 SðÞVr

;

where V is the volume of the reactor, F is the flowrate, S

is the substrate concentration in the feed

stream, and r

is the consumption rate of substrate. Growth of the microbial population X that feeds on

the substrate S is described by the kinetic rate equation

¼ μX;

where X is the cell mass, and the dependency of μ, the specific growth rate, on S is prescribed by the

Monod equation

μ ¼

max

þ S

;

where μ

max

is the maximum specific growth rate for X when S is not limiting. The substrate is consumed

at a rate equal to

XjS

S;max

þ S

;

where Y

X|S

is the called the yield factor and is defined as follows:

XjS

increase in cell mass

quantity of substrate consumed

415

7.1 Introduction

7.2 Euler’s methods

The simplest method devised to integrate a ﬁrst-order ODE is the explicit Euler

method. Other ODE integration methods may be more complex, but their underlying

theory parallels the algorithm of Euler’s method. Although Euler’s method is seldom

used, its algorithm is the most basic of all methods and therefore serves as a suitable

introduction to the numerical integration of ODEs. Prior to performing numerical

integration, the ﬁrst-order ODE should be expressed in the following format:

¼ ft; yðÞ; yt¼ 0ðÞ¼y

: (7:14)

ft; yðÞis the instantaneous slope of the function ytðÞat the point t ; yðÞ. The slope of y is,

in general, a function of t and y and varies continuously across the spectrum of values

and a < y < b. When no constraint is speciﬁed, inﬁnite solutions exist to the

ﬁrst-order ODE, dy=dt ¼ ft; yðÞ. By specifying an initial condition for y,i.e.

yt¼ 0ðÞ¼y

, one ﬁxes the initial point on the path. From this starting point, a

unique path can be traced for y starting from t = 0 by calculating the slope, ft; yðÞ,

at each point of t for 0

. The solution is said to be unique because, at any point t,

the corresponding value of y cannot have two simultaneous values. There are addi-

tional mathematical constraints that must be met to guarantee a unique solution for a

ﬁrst-order ODE, but these theorems will not be discussed here. We will assume in this

chapter that the mathematical requirements for uniqueness of a solution are satisﬁed.

A ﬁrst-order ODE will always be solved as an initial-value problem (IVP) because

only one constraint must be imposed to obtain a unique solution, and the

The predictive Tsuchiya equations of E. coli consumption by the amoeba Dictyostelium discoidium are

listed below (Bailey and Ollis, 1986; Tsuchiya et al., 1972). The limiting substrate is glucose. It is

assumed that the feed stream is sterile.

 SðÞ

;max

þ S

; (7:11)

¼

;max

þ S



;max

þ N

; (7:12)

¼

;max

þ N

: (7:13)

The kinetic parameters that define the predator–prey system are (Tsuchiya et al., 1972):

¼ 0:0625/hr;

¼ 0:5 mg/ml;

;max

¼ 0:25/hr, μ

;max

¼ 0:24/hr;

¼ 5  10

4

mg glucose/ml, K

¼ 4  10

bacteria/ml;

1=Y

¼ 3:3  10

10

mg glucose/bacterium;

1=Y

¼ 1:4  10

bacteria/amoeba.

At t =0, N

=13× 10

bacteria/ml and N

=4× 10

amoeba/ml. Integrate this system of

equations to observe the population dynamics of this predator–prey system. What happens to the

dynamics when a step change in the feed substrate concentration is introduced into the flow system?

416

Numerical integration of ODEs