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

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

that the transport characteristics of the drug (drug dissolution rate, diffusion of drug in three-

dimensional aqueous channels within the degrading polymeric matrix, diffusion of drug in the hydro-

dynamic layer surrounding the particle, polymer degradation kinetics, and kinetics of cleavage of

covalently attached drug from polymer) can be conveniently represented by a time-invariant mass

transfer coefficient k.

A mass balance on a drug-loaded PLGA microparticle when placed in an aqueous medium is given

by the first-order ODE



¼ kA C

 C

ðÞ;

where M is the mass of the drug remaining in the particle, k is the mass transfer coefficient that

characterizes the flux of solute per unit surface area per unit concentration gradient, A is the

instantaneous surface area of the sphere, C

is the solubility of the drug in water, and C

is the

concentration of drug in the medium. If the particle is settling under gravity in the quiescent medium,

the sedimentation rate of the particle is assumed to be slow. The concentration of drug immediately

surrounding the particle is assumed to be equal to the maximum solubility of drug at body

temperature.

If C

drug

is the uniform concentration of drug within the microsphere that does not change with time

during the polymer degradation process, and a is the instantaneous particle radius, we have

M ¼ C

drug

πa



The mass balance equation can be rearranged to obtain

¼

drug

 C

ðÞ: (7:62)

The size a of the microsphere shrinks with time.

For a spherical particle settling in Stokes flow under the effect of gravity, the momentum or force

balance that yields a first-order ODE in velocity was introduced in Example 7.2. However, the changing

particle size produces another term in the equation (adapted from Pozrikidis (2008)). The original

momentum balance for a settling sphere in Stokes flow,

πρ

dða

vÞ

þ 6μπav 

πa

ðρ

 ρ

Þg ¼ 0;

now becomes

þ 3a

¼ a

ðρ

 ρ

g 

9μa

2ρ

v: (7:63)

Substituting Equation 7.62 into Equation 7.63, we get

ðρ

 ρ

g 

9μ

v þ

3vk

drug

 C

ðÞ: (7:64)

The initial condition of Equation (7.64) is vð0Þ¼



t¼0

¼ 0. Equation (7.64) is nonlinear in a.

We are asked to determine the size, a

, of the largest drug-loaded microparticle that, when added to

the top of a chamber, will dissolve by the end of its travel down some vertical distance through a body of

fluid. Figure 7.15 illustrates the geometry of the system. Under gravity, the drug-loaded particle sinks in

low Reynolds number flow. As the particle falls through the medium, the drug is released in a controlled

and continuous manner.

To calculate the size of the particle that, starting from rest, will shrink to zero radius after falling a

distance L, we simultaneously integrate Equations (7.62) and (7.64). However, the time interval of

integration is unknown. Integration of the system of ODEs must include the first-order equation

467

7.7 Shooting method for boundary-value problems

¼ v

to calculate the distance traveled, so we can stop the integration when z ¼ L. It is more convenient to

integrate for a known interval length. Therefore, we change the independent variable of Equations (7.62)

and (7.64)fromt to z using the rule

Equations (7.62) and (7.64) become

¼

drug

 C

ðÞ; 0  z  L; (7:65)

ðρ

 ρ



9μ

drug

 C

ðÞ; 0  z  L: (7:66)

The boundary conditions for Equations (7.65) and (7.66) are

v 0ðÞ¼0; aLðÞ¼0:

We use the secant method to solve for a

. We try different values of a

as suggested by Equation (7.61)

to satisfy the boundary condition aLðÞ¼0.

The drug loading in a single microparticle is 20% w/w. The initial distribution of drug is uniform

within the particle. The parameters are assigned the following values at T ¼ 37 °C:

¼ 0:1 mg/ml (Sugano, 2008),

¼ 0:0 mg/ml,

drug

¼ 0:256 g/cm

k ¼ 1  10

6

cm/s (Siepmann et al., 2005),

¼ 1:28 g/cm

(Vauthier et al., 1999),

¼ 1:0 g/cm

Figure 7.15

Schematic diagram of the process of controlled drug release from a microparticle settling in aqueous fluid.

Solid core

Water

Air

Degraded hollow shell

468

Numerical integration of ODEs

μ ¼ 0:7cP,

L ¼ 10 cm.

The CGS system of units is used. However, the equations will be designed such that a

can be entered in

units of μ m. We make the following assumptions.

(1) As the polymer matrix-loaded drug dissolves into the aqueous surroundings, concomitantly the

outer particle shell devoid of drug completely disintegrates, or is, at the least, sufficiently hollow

that for all practical purposes it does not exert a hydrodynamic drag. The radius a therefore refers to

the particle core that has a density of ρ

¼ 1:28 g/cm

(2) The density, ρ

, of the particle core that contains the bulk of the drug does not change. In other

words, hydrolysis of polymer and dissolution of drug occur primarily at the shell surface and

penetrate inwards as a wavefront.

(3) The velocity is small enough that mass transfer is dominated by diffusional processes, and the

mass transfer coefficient is independent of particle velocity (Sugano, 2008).

On making the following substitutions

¼ a; y

¼ v; x ¼ z

into Equations (7.65) and (7.66), we obtain

¼

drug

 C

ðÞ; (7:67)

ðρ

 ρ



9μ

drug

 C

ðÞ: (7:68)

We will not be able to solve Equations (7.67) and (7.68) using the boundary conditions stated above.

Can you guess why? The solution becomes singular (blows up) at a = 0 and at v = 0. We must thus

make the boundary conditions more realistic.

(1) If we want 90% of the drug contained in the microsphere to be released, i.e. 90% of the sphere

volume to vanish during particle descent in the fluid, the boundary condition at x ¼ L becomes

x ¼ LðÞ¼0:464a

In other words, when the radius of the particle drops to below 46.4% of its original size, at least 90%

of the drug bound in the polymer matrix will have been released.

(2) We assume that the microparticle is introduced into the fluid volume from a height such that, when

the particle enters the aqueous medium, it has attained terminal Stokes velocity in air at T =20°C

(in air) ~ O 10

5



). Accordingly, the initial velocity at x = 0 is 0:016a

cm/s, where a

measured microns.

We also want to know the time taken for drug release. Let y

¼ t. We add a third equation to the ODE

set,

(7:69)

The secant equation (Equation (7.61)) is rewritten for this problem as

0;iþ1

¼ a

0;i



0;i

 a

0;i1



1;i

LðÞ0:464a

0;i



1;i

ðÞ

 0:464a

0;i



 y

1;i1

ðÞ

 0:464a

0;i1



(7:70)

Example 7.2 demonstrated that the problem of time-varying settling velocity is stiff. We use ode15s

to solve the IVP defined by Equations (7.67) – (7.69) and the initial conditions

469

7.7 Shooting method for boundary-value problems

x ¼ 0ðÞ¼a

; y

x ¼ 0ðÞ¼0:016a

; y

x ¼ 0ðÞ¼0:

Care must be taken to start with two guessed values for a

that are close to the correct initial condition

that satisfies the boundary conditions at x ¼ L . We choose the two guessed values for the initial

particle size to be a

0;1

¼ 1:2 μm and a

0;2

¼ 1:3 μm. MATLAB programs 7.16 and 7.17 contain the code

used to obtain the result.

After seven iterations of Equation (7.70), we obtain a solution that satisfies our tolerance condition:

¼ 1:14 μm. The time taken for 90% of the particle volume to dissolve is 43.6 hours. Figures 7.16 and

7.17 depict the change in radius and velocity of the particle as it traverses downward. Note that initial

velocity retardation is immediate and occurs over a very short time of O 10

6



s. This region of

integration is very stiff.

Figure 7.16

Change in microparticle size with distance traversed as the drug–polymer matrix erodes and dissolves into the

medium.

0 2 4 6 8 10

0.4

0.6

0.8

1.2

1.4

x (cm)

a (μm)

Figure 7.17

Change in microparticle settling velocity with distance traversed as the drug–polymer matrix erodes and dissolves

into the medium. (a) Large retardation in velocity of the particle over the first 0.00013 μm of distance traveled.

(b) Gradual velocity decrease taking place as the particle falls from 0.00013 μmto10μm.

(a) (b)

–8

2 × 10

−8

100

150

200

250

x (cm)

v (μm/s)

1.3e−008 5 10

0.2

0.4

0.6

0.8

1.2

x (cm)

v (μm/s)

470

Numerical integration of ODEs

MATLAB program 7.16

% Use the shooting method to solve the controlled drug delivery

problem.

clear all

% Variables

maxloops = 50; % maximum number of iterations allowed

L = 10; % length of interval (cm)

y0(3) = 0; % initial time (s)

% First guess for initial condition

a0_1 = 1.2; % radius of particle (um)

y0(1)= a0_1;

y0(2) = 0.016*a0_1^2; % velocity (cm/s)

% Call ODE Solver: First guess

[x, y] = ode15s(‘drugdelivery’, L, y0);

y1(1) = y(end,1); % End-point value for y1

% Second guess for initial condition

a0_2 = 1.3; % radius of particle (um)

y0(1) = a0_2;

y0(2) = 0.016*a0_2^2; % velocity (cm/s)

% Call ODE Solver: Second guess

[x, y] = ode15s(‘drugdelivery’, L, y0);

y1(2) = y(end,1); % End-point value for y1

% Use secant method to ﬁnd correct initial condition

for i = 1:maxloops

% Secant Equation

a0_3 = a0_2 - (a0_2 - a0_1)*(y1(2)-a0_2*0.464)/ ...

((y1(2) - a0_2*0.464) - (y1(1) - a0_1*0.464));

% New guess for initial condition

y0(1) = a0_3;

y0(2) = 0.016*a0_3^2; % velocity (cm/s)

% Solving system of ODEs for improved value of initial condition

[x, y] = ode15s(‘drugdelivery’, L, y0);

y1(3) = y(end,1);

% Check if ﬁnal value satisﬁes boundary condition

if abs(y1(3) - a0_3*0.464) < a0_3*0.001 % tolerance criteria

break

end

a0_1 = a0_2;

a0_2 = a0_3;

y1(1) = y1(2);

y1(2) = y1(3);

end

MATLAB program 7.17

function f = drugdelivery(x, y)

% Evaluate f(t,y) for controlled-release drug delivery problem

% Constants

471

7.7 Shooting method for boundary-value problems

7.8 End of Chapter 7: key points to consider

(1) An ordinary differential equation (ODE) contains derivatives with respect to only one

independent variable.

(2) An ODE of order n requires n constraints or boundary values to be speciﬁed before a

unique solution can be determined. An ODE is classiﬁed based on where the con-

straints are speciﬁed. An initial-value problem has all n constraints speciﬁed at the

starting point of the interval and a boundary-value problem has some conditions

speciﬁed at the beginning of the interval and the rest speciﬁed at the end of the interval.

(3) To integrate an ODE numerically, one must divide the integration interval into

many subintervals. The number of subintervals depends on the step size and the size

of the integration interval. The endpoints of the subintervals are called nodes or time

points. The solution to the ODE or system of ODEs is determined at the nodes.

(4) A numerical integration method that determines the solution at the next time point,

kþ1

, using slopes and/or solutions calculated at the current time point t

, and/or

prior time points (e.g. t

k1

; t

k2

) is called an explicit method, e.g. Euler’s forward

method. A numerical integration method whose difference equation requires a slope

evaluation at t

kþ1

in order to ﬁnd the solution at the next time point, t

kþ1

, is called an

implicit method, e.g. Euler’s back ward method. The unknown terms in an implicit

formula are located on both sides of the difference equation.

(5) The local truncation error at any time point in the integration interval is the numerical

error due to the truncation of terms, if one assumes that the previous solution(s) used to

calculate the solution at the next time point is exact. The global truncation error at any

time point in the integration interval is the error between the numerical approximation

and the exact solution. The global error is the sum of (i) accumulated errors from

previous time steps, and (ii) the local truncation error at the current integration step.

(6) The order of a numerical ODE integration method is determined by the order of the

global truncation error. An nth-order method has a global truncation error of Oh

ðÞ.

(7) A one-s tep integration method uses the solution at the current time point to predict

the solution at the next time point. A multistep integration method evaluates the slope

function at previous time steps and collectively uses the knowledge of slope behavior

at multiple time points to determine a solution at t

kþ1

(8) An ODE or system of ODEs is said to be stable if its solution is bounded at all time t.

This is possible if the real parts of all eigenvalues of the Jacobian matrix are negative.

If even one eigenvalue has a positive real part, the solut ion will grow in time without

Cs = 0.1e-3; % solubility of drug (g/m1)

Cm = 0.0; % concentration of drug in medium (g/m1)

Cdrug = 0.256; % concentration of drug in particle (g/cm^3)

rhop = 1.28; % density of particle (g/cm^3)

rhof = 1.0; % density of medium (g/cm^3)

mu = 0.7e-2; % viscosity of medium (Poise)

k = 1e-6; % mass transfer coefﬁcient (cm/s)

g = 981; % acceleration due to gravity (cm/s^2)

f = [-k/y(2)/Cdrug*(Cs - Cm)*1e4; ...

(rhop-rhof)/rhop*g/y(2) - 9*mu/2/(y(1)*1e-4)^2/rhop + ...

3*k/(y(1)*1e-4)/Cdrug*(Cs - Cm); ...

1/y(2)];

472

Numerical integration of ODEs

bound. A numerical solution is said to be stable when the error remains constant or

decays with time, and is unstable if the error grows exponentially with time. The

stability of a numerical solution is dependent on the na ture of the ODE system, the

properties of the ODE solver, and the step size chosen. Implicit methods have better

stability properties than explicit methods.

(9) A stiff differential equation is difﬁcult to integrate since it has one or more rapidly

vanishing or accelerating transients contained in the solution. To integrate a stiff

ODE, one must use a numerical scheme that has superior stability properties. An

integration method with poor stability properties must necessarily use vanishingly

small step sizes to ensure that the numerical solution remains stable. When the step

size is too small, round-off error will preclude a stable solution.

(10) The shooting method is used to solve boundary-value ODE problems. This method

converts the boundary-value problem into an initial-val ue problem by guessing

values for the unspeciﬁed initial conditions. The numerical estimates for the boun-

dary endpoint values at the end of the integration are compared to the target values.

One can use a nonlinear root-ﬁnding algorithm to improve iterations of the guessed

initial condition(s).

(11) Table 7.5 lists the MATLAB ODE functions that are discussed in this chapter.

7.9 Problems

7.1. Monovalent binding of ligand to cell surface receptors Consider a monovalent

ligand L binding reversibly to a monovalent receptor R to form a receptor/ligand

complex C:

R þ L

!



where C represents the concentration of receptor (R) – ligand (L) complex, and k

and k

are the forward and reverse reaction rates, respectively. If we assume that the

Table 7.5. MATLAB ODE functions discussed in the chapter

ODE

function Algorithm Type of problem

ode45 RKF method: RK−4 and RK−5 pair

using the Dormand and Prince

algorithm

non-stiff ODEs

ode23 RKF method: RK−2 and RK−3 pair

using the Bogacki and Shampine

algorithm

non-stiff ODEs when high

accuracy is not required

ode113 Adams–Bashforth–Moulton method non-stiff ODEs when high

accuracy is required

ode15s multistep method stiff problems; use on the ﬁrst try

ode23s one-step method stiff problems at low accuracy

ode23t trapezoidal rule stiff equations

ode23tb implicit RK formula stiff equations at low accuracy

473

7.9 Problems

ligand co ncentration remains constant at L

, the following differential equation

governs the dynamics of C (Lauffenburger and Linderman, 1993):

¼ k

 C½L

 k

where R

is the total number of receptor molecules on the cell surface.

(a) Determine the stability criterion of this differential equation by calculating the

Jacobian. Supposing that we are interested in studying the binding of ﬁbronectin

receptor to ﬁbronectin on ﬁbroblast cells (R

=5× 10

sites/cell, k

=7× 10

(M min), k

= 0.6/min). What is the range of ligand concentrations L

for which

this differential equation is stable?

(b) Set up the numerical integration of the differential equation using an impl icit

second-order rule, i.e. where the error term scales as h

and h = 0.1 min is the

step size. Use L

¼ 1 μM. Perform the integration until equilibrium is reached.

How long does it take to reach equilibrium (99% of the ﬁnal concentration C is

achieved)? (You will need to calculate the equilibrium concentration to perform

this check.)

7.2. ODE model of convection–diffusion–reaction The following equation is a differ-

ential model for convection, diffusion, and reaction in a tubular reactor, assuming

ﬁrst-order, reversible reaction kinetics:



¼ Da  C;

where Pe is the Peclet number, which measures the relative importance of convection

to diffusion, and Da is the Damkohl er number, which compares the reaction time

scale with the characteristic time for convective transport.

Perform a stability analysis for this equation.

(a) Convert the second-order system of equations into a set of coupled ﬁrst-order

ODEs. Calculate the Jacobian J for this system.

(b) Find the eigenvalues of J. (Find the values of λ such that A  λI

¼ 0.)

7.3. Chlorine loading of the stratosphere Chloroﬂuorocarbons (CFCs), hydroﬂuor-

ocarbons (HCFCs), and other haloalkanes are well known for their ozone-depleting

characteristics. Accordingly, their widespread use in industries for refrigeration, ﬁre

extinguishing, and solvents has been gradually phased out. CFCs have a long life-

span in the atmosphere because of their relative inertness. These compounds enter

the stratosphere, where UV radiation splits the molecules to produce free chlorine

radicals (Cl·) that catalyze the destruction of ozone (O

) located within the strato-

sphere. As a result, holes in the ozone have formed that allow harmful UV rays to

reach the Earth’s surface.

CFC molecules released into the atmosphere cycle every three years or so between

the stratosphere and the troposphere. Halogen radicals are almost exclusively pro-

duced in the stratosphere. Once transp orted back to the troposphere, these radicals

are usually washed out by rain and return to the soil (see Figure P7.1).

The HCFCs are not as long-lived as CFCs since they are more reactive; they break

down in the troposphere and are eliminated after a few years. The concentrations of

halocarbons in the stratosphere and troposphere as well as chlorine loading in the

stratosphere can be modeled using a set of ordinary differential equations (Ko et al.,

1994). If B

is the concentration of undissociated halocarbons in the stratosphere, B

is the concentration of undissociated halocarbons in the troposphere, and C is the

474

Numerical integration of ODEs

quantity of chlorine in the stratosphere, then we can formulate the following

equations:

 fB



;

 B



;



where

τ is the characteristic time for replacing stratospheric air with tropospheric air,

f = 0.15/0.85 is the scaled fraction of tropospher ic air that enters the stratosphere

(15% of atmospheric mass lies in the stratosphere),

is the halocarbon lifetime in the troposphere that characterizes the time until

chemical degradation,

is the halocarbon lifetime in the stratosphere that characterize s the time until

chemical degradation.

It is assumed that chlorine is quickly lost from the troposphere. The parameter

values are set as τ ¼ 3 yr, L

= 1000 yr, and L

= 5 yr for the halocar bon CFCl

(Ko

et al., 1994).

(a) What is the Jacobian for this ODE system? Use the MATLAB function eig to

calculate the eigenvalues for J. Study the magnitude of the three eigenvalues,

and decide if the system is stiff. Accordingly, choo se the appropriate ODE

solver.

(b) If 100 kg of CFC is released at time t = 0, what is the CFC loading in the

troposphere and chlorine loading within the stratosphere after 10 years? And

after 100 years?

7.4. Population dynamics among hydrocarbon-consuming bacterial species Pseudomonas

is a class of bacteria that converts methane in the gas phase into methanol.

Secreted methanol participates in a feedback inhibition control process inhibiting

further growth. When grown with another bacterial species, Hyphomicrobium , both

species beneﬁt from an advantageous interdependent interaction (mutualism).

Hyphomicrobium utilizes methanol as a substrate and thereby lowers the methanol

concentration in the surrounding medium promoting growth of Pseudomonas.

Wilkinson et al.(1974) investigated the dynamic behavior of bacterial growth of

these two species in a chemostat. Their developments are reproduced here based on

the discussion of Bailey and Ollis (1986). Let x

denote the Pseudomonas population

mass and let x

denote the Hyphomicrobium population mass. Dissolved oxygen is a

Figure P7.1

First two atmospheric layers.

Troposphere

Stratosphere

2Cl·

Ozone layer

Sea level

475

7.9 Problems

requirement for Pseudomonas growth, and its growth rate can be described by the

following kinetic rate equation:

;max

þ c



1 þ S=K

This rate equation takes into account the inhibitory action of methanol (S), where K

is the inhibition constant. The growth of Pseudomonas exhibits a monod dependence

on dissolved oxygen concentration, c

, and S is the methanol concentration.

The growth of Hyphomicrobium is not oxygen-limited, but is substrate (methanol)-

dependent, and its growth kinetics are described by the following rate equation:

;max

þ S

The reactions take place in a chemostat (see Figure P7.2, and see Box 7.4A for an

explanation of this term).

The following mass balance equations describe the time-dependent concentra-

tions of the substrates and two bacterial populations:

¼ k

 c

ðÞ



;

where k

 c

ðÞis the mass transfer rate of oxygen at the gas–liquid interface.

We have

¼

S þ



;

¼

þ r

;

¼

þ r

The reaction parameters estimated by Wilkinson et al.(1974) are listed below.

;max

¼ 0:185/hr; μ

;max

¼ 0:185/hr,

¼ 1  10

5

g/l; K

¼ 5  10

6

g/l,

¼ 1  10

4

g/l,

¼ 0:2gpseudomonas bacteria/g O

consumed,

¼ 5:0gpseudomonas bacteria produced/g methanol produced,

Figure P7.2

Flow into and out of a chemostat.

Stirrer

476

Numerical integration of ODEs