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

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

2 Systems of linear equations

2.1 Introduction

A linear equation is one that has a simple linear dependency on all of the variables or

unknowns in the equation. A generalized linear equation can be written as:

þ a

þþa

¼ b:

The variables of the equation are x

; x

; ...; x

. The coefﬁcients of the variables

; a

; ...; a

, and the constant b on the right-hand side of the equation represents the

known quantities of the equation. A system of linear equations consists of a set of two

or more linear equations. For example, the system of equations

þ x

¼ 4;

þ 3x

 4x

¼ 4

forms a set of two linear equ ations in three unknowns. A solution of a system of

linear equations consists of the values of the unknowns that simultaneously satisfy

each linear equation in the system of equations. As a rule of thumb, one requires n

linear equations to solve a system of equations in n unknowns to obtain a unique

solution. A problem is said to have a unique solution if only one solution exists that

can exactly satisfy the governing equations. In order for a system of linear equations

in three unknowns to yield a unique solution, we must have at least three linea r

equations. Therefore, the above set of two equations does not have a unique

solution. If we complete the equation set by including one more equation, then we

have three equations in three unknowns:

þ x

¼ 4;

þ 3x

 4x

¼ 4;

 3x

¼ 0:

The solution to this system of equations is unique: x

¼ 1; x

¼ 2; x

¼ 1.

Sometimes one may come across a situation in which a system of n equations in n

unknowns does not have a unique solution, i.e. more than one solution may exist

that completely satisﬁes all equations in the set. In other cases, the linear problem

can have no solution at all. Certain properties of the coefﬁcients and constants of the

equations govern the type of solution that is obtained.

Many standard engineering problems, such as force balances on a beam or truss,

ﬂow of a reagent or chemical through a series of stirred tanks, or ﬂow of electricity

through a circuit of wires, can be modeled using a system of linear equations.

Solution methods for simultaneous linear equations (that contain two or more

unknowns) are in most cases straightforward, intuitive, and usually easier to

implement than those for nonlinear systems of equations. Linear equations are

solved using mathematical techniques derived from the foundations of linear algebra

theory. Mastery of key concepts of linear algebra is important to apply solution

methods developed for systems of linear equations, and understand their limitations.

Box 2.1A Drug development and toxicity studies: animal-on-a-chip

Efficient and reliable drug screening methods are critical for successful drug development in a cost-

effective manner. The FDA (Food and Drug Administration) mandates that the testing of all proposed

drugs follow a set of protocols involving (expensive) in vivo animal studies and human clinical trials. If

the trials are deemed successful, only then may the approved drug be sold to the public. It is most

advantageous for pharmaceutical companies to adopt in vitro methods that pre-screen all drug

candidates in order to select for only those that have the best chances of success in animal and

human trials. Key information sought in these pre-clinical trial experiments include:

absorption rate into blood and tissues, and extent of distribution of drug in the body,

metabolism of the drug,

toxicity effects of the drug and drug metabolites on the body tissues, and

routes and rate of excretion of the drug.

The biochemical processing o f drug by the body can be represented by physiologically based

mathematical models that describe the body’s response to the drug. These models are “physio-

logical” because the parameters of the equations are obtained from in vivo animal experiments and

therefore reliably predict the physiological response of the body when the d rug is introd uced.

Mathematical models that predict the distribution (concentration of drug in various par ts of the

body), metabolism (breakdown of drug by enzymatic reactions), covalent binding of drug o r drug

metabolites with ce llular components or extracellular proteins, and elimination of the drug from the

body are called pharmacokinetic models.SeeBox 2.3 for a discussion on “pharma cok ine tic s” as a

branch of pharmacology.

A physiologically based pharmacokinetic model (PBPK) using experimental data from in vitro

cell-based assays was successfully developed by Quick and Shuler (1999) at Cornell University. This

model could predict many of the drug (toxicant: naphthalene) responses observed in vivo, i.e.

distribution, enzymatic consumption of drug, and non-enzymatic binding of drug metabolites to tissue

proteins, in several different organs of the rat and mouse body (e.g. lung, liver, and fatty tissue). The

merit in developing in vitro-based physiological models stems from the fact that data from in vitro

studies on animal cells are relatively abundant; in vivo animal studies are considerably more expensive

and are subject to ethical issues. Performing experimental studies of drug effects on humans is even

more difficult. Instead, if an overall human body response to a drug can be reliably predicted by

formulating a PBPK model based on samples of human cells, then this represents advancement by

leaps and bounds in medical technology.

A novel microfabricated device called microscale cell culture analog (μCCA) or ‘animal-on-a-

chip’ was recently invented by Shuler and co-workers (Sin et al., 2004). The μCCA uses microchannels

to transport media to different compartments that mimic body organs and contain representative animal

cells (hence “animal-on-a-chip”). The flow channels and individual organ chambers are etched onto a

silicon wafer using standard lithographic methods. The μCCA serves as a physical replica of the PBPK

model and approximates the distribution, metabolism, toxic effects, and excretion of drug occurring

within the body. The advantages of using a μCCA as opposed to a macro CCA apparatus (consisting of a

set of flasks containing cultured mammalian cells, tubes, and a pump connected by tubing) are several-

fold:

(1) physiological compartment residence times, physiological flowrates (fluid shear stresses), and

liquid-to-cell ratios are more easily realized;

Systems of linear equations

(2) rate processes occur at biological length scales (on the order of 10 μm), which is established only

by the microscale device (Sin et al., 2004);

(3) smaller volumes require fewer reagents, chemicals and cells (especially useful when materials are

scarce); and

(4) the device is easily adaptable to large-scale automated drug screening applications.

Figure 2.1 is a simplified block diagram of the μCCA (Sin et al., 2004) that has three compartments:

lung, liver, and other tissues.

The effects of the environmental toxicant naphthalene (C

) – a volatile polycyclic aromatic

hydrocarbon – on body tissues was studied by Shuler and co-workers using the μCCA apparatus

(Viravaidya and Shuler, 2004; Viravaidya et al., 2004). Naphthalene is the primary component of moth

balls (a form of pesticide), giving it a distinctive odor. This hydrocarbon is widely used as an

intermediate in the production of various chemicals such as synthetic dyes. Exposure to large amounts

of naphthalene is toxic to humans, and can result in conditions such as hemolytic anemia (destruction

of red blood cells), lung damage, and cataracts. Naphthalene is also believed to have carcinogenic

effects on the body.

On exposure to naphthalene, cytochrome P450 monooxygenases in liver cells and Clara lung cells

catalyze the conversion of naphthalene into naphthalene epoxide (naphthalene oxide). The circulat-

ing epoxides can bind to glutathione (GSH) in cells to form glutathione conjugates. GSH is an

antioxidant and protects the cells from free radicals, oxygen ions, and peroxides, collectively known as

reactive oxygen species (ROS), at the time of cellular oxidative stress. Also, GSH binds toxic

metabolic by-products, as part of the cellular detoxification process, and is then eliminated from the

body in its bound state (e.g. as bile). Loss of GSH functionality can lead to cell death.

Naphthalene epoxide is consumed in several other pathways: binding to tissue proteins, sponta-

neous rearrangement to naphthol, and conversion to naphthalene dihydrodiol by the enzyme

epoxide hydrolase. Both naphthol and dihydrodiol are ultimately converted to naphthalenediol,

which then participates in a reversible redox reaction to produce naphthoquinone. ROS generation

occurs in this redox pathway. The build up of ROS can result in oxidative stress of the cell and cellular

damage. Also, naphthoquinone binds GSH, as well as other proteins and DNA. The loss and/or alteration

Figure 2.1

Block diagram of the μCCA (Sin et al., 2004). The pump that recirculates the media is located external to the silicon

chip. The chip houses the organ compartments and the microchannels that transport media containing drug

metabolites to individual compartments.

Lung

Other tissues

Liver

Out

Pump

Silicon chip

2.1 Introduction

of protein and DNA function can result in widespread havoc in the cell. Figure 2.2 illustrates the series of

chemical reactions that are a part of the metabolic processing of naphthalene into its metabolites.

Naphthalene is excreted from the body in the form of its metabolites: mainly as GSH conjugates,

naphthalene dihydrodiol, and naphthol. Quick and Shuler (1999) developed a comprehensive PBPK

model for exposure to naphthalene in mice. In mice, exposure to naphthalene results in cell death in the

lung.

Here, we envision the steady state operation of Shuler’s μCCA that might mimic continuous exposure

of murine cells to naphthalene over a period of time. Figure 2.3 is a material balance diagram for

naphthalene introduced into the μCCA depicted in Figure 2.1. The feed stream contains a constant

concentration of naphthalene in dissolved form. Conversion of naphthalene to epoxide occurs in the

lung and liver in accordance with Michaelis–Menten kinetics (see Box 4.3 for a discussion on

Michaelis–Menten kinetics). Both compartments are treated as well-mixed reactors. Therefore, the

prevailing concentration of naphthalene within these compartments is equal to the exit concentration.

No conversion of naphthalene occurs in the other tissues (ot). Naphthalene is not absorbed anywhere in

the three compartments. The purge stream mimics the loss of the drug/toxicant from the circuit via

excretory channels or leakage into other body parts (see Figure 2.3).

In engineering, a mass balance is routinely performed over every compartment (unit) to account for

transport, distribution among immiscible phases, consumption, and generation of a particular product,

drug, or chemical. A mass balance is a simple accounting process of material quantities when

phenomena such as convective transport and/or physical processes (diffusion, absorption) and/or

chemical processes (reaction, binding) are involved. A generic mass or material balance of any

component A over a unit, equipment, or compartment is represented as:

ðflow in of A  flow out of AÞþðgeneration of A  consumption of AÞ

¼ accumulation of A in unit:

Figure 2.2

Naphthalene metabolic pathway in mice (Viravaidya et al., 2004).

Naphthalene

epoxide

Naphthalene

dihydrodiol

Naphthol Epoxide–GSH

conjugates

Cytochrome P450

Epoxide

hydrolase

Molecular

rearrangement

Glutathione S-transferase

Naphthalenediol

Cytochrome

P450

Enzymatic

conversion

Naphthoquinone

ROS generation

Naphthoquinone–GSH

conju

ates

Systems of linear equations

In steady state processes, the concentrations of components as well as the volumetric flowrates

do not change with time, and thus there is no accumulation of material anywhere in the system. For a

steady state process, the material balance is:

ðflow in of A  flow out of AÞþðgeneration of A  consumption of AÞ¼0:

Mass balances over individual compartments (units) as well as an overall mass balance over the

entire system of units provide a set of equations that are solved to obtain the unknowns such as the

outlet concentrations of the species of interest A from each unit.

On performing a mass balance of naphthalene over the two chambers – lung and liver – we obtain

two equations for the unknowns of the problem, C

lung

and C

liver

Lung compartment

feed

þ RQ

liver

þ Q

lung





max;P450-lung

lung

m;P450-lung

þ C

lung

 Q

lung

¼ 0: (2:1)

Liver compartment

liver

lung



max;P450-liver

liver

m;P450-liver

þ C

liver

 Q

lung

¼ 0(2:2)

The mass balance of naphthalene over the other tissues compartment is trivial since no loss or

conversion of mass occurs in that compartment.

The values for the parameters in Equations (2.1) and (2.2) are provided below.

lung

: flowrate through lung compartment = 2 μl/min;

liver

: flowrate through liver compartment = 0.5 μl/min;

Figure 2.3

Material balance of naphthalene in the μCCA.

Other tissues

feed

, C

feed

lung

, C

lung

, C

lung

liver

, C

lung

liver

, C

liver

, C

lung

liver

R(Q

liver

), C

out

Recycle stream

Additional equations:

Exit stream from lung is split into two streams: Q

lung

= Q

liver

+ Q

Flow into chip = Flow out of chip: Q

feed

= (1

–

R)Q

lung

Well-mixed lung

reactor

Well-mixed liver

reactor

Purge stream

(1-R)(Q

liver

), C

out

2.1 Introduction

: flowrate through other tissues compartment = 1.5 μl/min;

feed

: concentration of naphthalene in feed stream = 390 μM;

feed

: flowrate of feed stream = ð1  RÞQ

lung

;

max;P450-lung

: maximum reaction velocity for conversion of naphthalene into napthalene epoxide by

cytochrome P450 in lung cells = 8.75 μM/min;

m;P450-lung

: Michaelis constant = 2.1 μM;

max;P450-liver

: maximum reaction velocity for conversion of naphthalene into napthalene epoxide by

cytochrome P450 in liver cells = 118 μM/min;

m;P450-liver

: Michaelis constant = 7.0 μM;

lung

: volume of lung compartment = 2 mm  2mm 20 μm ¼ 8  10

8

l ¼ 0:08μ

liver

: volume of liver compar tment = 3:5mm 4:6mm 20 μm ¼ 0:322 μl

R: specifies the fraction of the exiting stream that re-enters the microcircuit.

Substituting these values into Equations (2.1) and (2.2), we obtain the following.

Lung compartment

780ð1  RÞþR 0:5C

liver

þ 1 :5C

lung





8:75C

lung

2:1 þ C

lung

0:08  2C

lung

¼ 0:

Liver compartment

0:5C

lung



118C

liver

7:0 þ C

liver

0:322  0:5C

liver

¼ 0:

Here, we make the assumption that C

lung

 K

m;P450-lung

and C

liver

 K

m;P450-liver

. As a result, the

Michaelis–Menten kinetics tends towards zeroth-order kinetics. We then have the following.

Lung compartment

780ð1  RÞþR 0:5C

liver

þ 1 :5C

lung



 0:7  2C

lung

¼ 0:

Liver compartment

0:5C

lung

 38  0:5C

liver

¼ 0;

which can be simplified further as follows.

Lung compartment

2  1:5RðÞC

lung

 0 :5RC

liver

¼ 779:3  780R: (2:3)

Liver compartment

lung

 C

liver

¼ 76: (2:4)

We have two simultaneous linear equations in two variables. Equations (2.3) and (2.4) will be solved

using solution techniques discussed in this chapter.

Some key assumptions made in setting up these equations are:

Systems of linear equations

2.2 Fundamentals of linear algebra

This section covers topics on linear algebra related to the solution methods for simulta-

neous linear equations. In addition to serving as a review of fundamental concepts of

linear algebra, this section introduces many useful MATLAB operators and functions

relevant to matrix construction and manipulation. Applicati on of mathematical tools

from linear algebra is common in many ﬁelds of engineering and science.

MATLAB or “MATrix LABoratory” is an ideal tool for studying and perform-

ing operations involving matrix algebra. The MATLAB computing environment

stores e ach data set in a matrix or multidimensional array format. For example, a

scalar value of a=1is stored as a 1 × 1 matrix. A character set such as ‘computer’

is treated as a 1 × 8 matrix of eight character elements. Dimensioning of arrays prior

to use is not required by MATLAB. Any algebraic, relational, or logical operation

performed on any data or variables is treated as a matrix or array operation.

MATLAB contains many built-in tools to analyze and manipulate matrices.

2.2.1 Vectors and matrices

In linear algebra, a vector is written as a set of numbers arranged in sequential

fashion either row-wise or column-wise and enclosed within brackets. The individual

entries contained in a vector are the elements. For example, a vector v that contains

three elements 1, 2 and 3 is written as:

v ¼ 123½or v ¼

1  3 matrix 3  1 matrix

It is not necessary to include a comma to delineate each entry within the brackets.

The vector v written row-wise is called a row vector and can be visualized as a matrix

of size 1 × 3, where 1 is the number of rows and 3 is the number of columns. On the

other hand, v is a column vector when its elements are arranged in a single column

and can be regarded as a matrix of size 3 × 1, where 3 is the number of rows and 1 is

the number of columns. An n-element column vector v is represented as follows:

v ¼

;

where v

are the individual elements or components of v and i =1, 2, ..., n. Each

element has a single subscript that speciﬁes its position in the sequence. Throughout

(1) all reactions are irreversible;

(2) all compartments are well-mixed;

(3) no density changes occur in the system;

(4) The tissue : medium partition coefficient is 1, i.e. there is an equal distribution between tissue and

medium.

2.2 Fundamentals of linear algebra

this text, variables representing vectors are denoted by lower-case bold roman

letters, and elements of a vector will be represented by subscripted lower-case

roman letters. Henceforth, every vector is assumed to be a column vector (unless

noted otherwise) as per common convention used by many standard linear algebra

textbooks.

A matrix is deﬁned as a two-dimensional rectangular array of numbers. Examples

of matrices are shown below:

A ¼

472

518



; B ¼



; C ¼

Matrix variables are denoted by bold upper-case roman letters in order to provide a

contrast with vectors or one-dimensional matrices. The size of a matrix is determined

by the num ber of rows and columns that the matrix contains. In the examples above,

A and C are rectangular matrices of sizes 2 × 3 and 3 × 2, respectively, and B is a

square matrix of size 2 × 2, in which the number of rows equal the number of

columns. An m × n matrix A, where m is the number of rows and n is the number

of columns, can be written in a generalized form at as shown below, in which every

element has two subscripts to specify its row and column position. The ﬁrst su bscript

speciﬁes the row location and the second subscript indicates the column position of

an element:

A ¼

... a

An elemen t in the ith row and jth column of A is represented as a

Using MATLAB

A vector or matrix can be created by using the matrix constructor operator, [ ].

Elements of the matrix are entered within square brackets and individual rows are

distinguished by using a semi-colon. Elements within a row can be separated using

either a comma or space. Here, we construct a 3 × 3 matrix A in MATLAB:

A=[123;456;789]

MATLAB provides the following output (since we have not suppressed output

generation by adding a semi-colon at the end of the statement):

123

456

789

There are several MATLAB functions available to create special matrices for

general use. The zeros function produces a matrix consisting entirely of zeros of

speciﬁed size. Similarly, the ones function produces a matrix consisting en tirely of

ones of speciﬁed size. For example,

Systems of linear equations

Z = zeros(2,3)

000

M = ones(4,4)

1111

Two other important matrix generating functions are the eye and diag func-

tions. The eye function produces the identity matrix I, in which the elements are

ones along the diagonal and zeros elsewhere. The diag function creates a diagonal

matrix from an input vector by arranging the vector elements along the diagonal.

Elements along the diagonal of a matrix can also be extracted using the diag

function. Exa mples of the usage of eye and diag are shown below.

I = eye(3)

100

010

001

D = diag([1 2 3])

100

020

003

b = diag(I)

To access a particular element of a matr ix, the row and column indices need to be

speciﬁed along with the matrix name as shown below. To access the element in the

second row and second column of matrix

123

456

789

we type

A(2,2)

ans =

If multiple elements need to be accessed simultaneously, the colon operator (:) is

used. The colon operator represents a series of numbers; for example, 1:n represents

a series of numbers from 1 to n in increments of 1:

2.2 Fundamentals of linear algebra

A(1:3, 2)

ans =

If access to an entire row or c olumn is desired, the colon by itself can be used as a

subscript. A stand-alone colon refers to all the elements along the speciﬁed column

or row, such as that shown below for column 2 of matrix A:

A(:,2)

ans =

MATLAB provides the size function that returns the length of each dimension

of a vector or matr ix:

[rows, columns] = size(A)

rows =

columns =

2.2.2 Matrix operations

Matrices (and vectors) can undergo the usual arithmetic operations of addition and

subtraction, and also scalar, vector, and matrix multiplication, subject to certain

rules that govern matrix operations. Let A be an m × n matrix and let u and v be

column vectors of length n × 1.

Addition and subtraction of matrices

Addition or subtraction of two matrices is carried out in element-wise fashion and

can only be performed on matrices of identical size. Only elements that have the

same subscripts or position within their respective matrices are added (or sub-

tracted). Thus,

A þ B ¼ a

þ b



mn

; where 1  i  m and 1  j  n:

Scalar multiplication

If a matrix A is multiplied by a scalar value or a single number p, all the elements of

the matrix A are multiplied by p:

pA ¼ pa



mn

; where 1  i  m and 1  j  n:

Using MATLAB

If we deﬁne

A = 3*eye(3)

Systems of linear equations