MSU

LIBRARIES

 

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES wi11
be charged if book is
returned after the date

stamped below.

 

 

 

 

NON-EQUILIBRIUI TEEENDDYNAIIC
IDDELING AND PARAMETER ESTIMATION OF
PEENOIENOLOGICAL COEFFICIENTS DESCRIBING COUPLED
TRANSPORT ACROSS A HEBRANE

By

Steven Patrick Nowlen

A.IEESIS

Submitted to
liehigen Stete University
in pertiel fullfillnent of the require-ents
for the degree of

ULSTER OF SCIENCE

Depertnent of Ieehenieel Engineering

1983

ABSTRACT

NON-EQUILIBRIUI'TEERIODTNAMIC
MODELING AND PARAIETER ESTIMKTION OF
PRENOIENOLOGICAL COEFFICIENTS DESCRIBING COUPLED
TRANSPORT ACROSS A MEMBRANE

By

Steven Patrick Nowlen

A diffusion chamber microscope stage has been developed which
subjects a email smmple, of cells (lOul) to a psuedo-step-change in
extracellular concentration of permeable and/or impermeable solutes.
The response of an individual cell to the induced osmotic imbalance
was documented via a series of photomicrographs. These images were
processed using simple image analysis techniques to yield the
volume-time history of the c311 of interest. These volume-time data
were used to estimate the values of the permeability parameters for
the cell membrane through one of three modified ordinary least squares
parameter estimation. methods linked to an irreversible thermodynamic
model of the membrane transport process. Data were successfully pro-
cessed for hamster embryos for both the single parameter omnotic
shrinkage. and the three parameter binary flow shrinkrswell cases of
the Iedmm and Ketchalsky permeability model. These results are con-

sistent with results reported previously by other investigators.

ACKNOWLEDGMENTS

I would like to thank my adviser. John McGrath. for his help and

guidance as well as for his friendship. I would also like to thank

Mohsen Shabana for providing much of the data presented as a part of

this work. and the other members of the RT? lab group for their help.
I would also like to express my appreciation to the Department of

Mechanical Engineering and to the Department of Engineering Research

whose finacial support made my graduate work possible.

TABLE OF CONTENTS

List of Tables

List of Figures

Nomenclature Listing

Chapter 1:
1.1:
1.2:

Chapter 2:
2.1:
2.2:
2.3:
2.4:
2.5:

2.6:

2.7
2.8:

Chapter 3:
3.1:
3.2:

Chapter 5:
5.1:
5.2:
Chapter 6:

Chapter 7:

Introduction
Motivation for Present Work
Statement of Objectives

Historical Deve10pment of Permeability Models
Historical Overview
Jacobs' Model
The K-K Permeability Formulation
The Kr! Resistance Formulation
Comparison of the Re! Resistance and Permeability
Formulations
A Power Series Solution to the K-K Permeability
Equations

: The Papanek Model

Nondimensionalization of the KéK Permeability
Equations

Parameter Estimation Methods

Parameter Estimation Overview

Ordinary Least Squares and Maximum Likelihood
Estimators

: The Gauss Minimization Method
: Box-Kanemasu Interpolation Method
: Marquardt's Method

: Data Gathering Techniques
: The Diffusion Chamber
: Characterization of the Dialysis Mombrane

Permeability

: Laboratory Tbchniques
: Image Processing Techniques

Data and Results
Osmotic Shrinkage of Unfertilized Hamster Ova
Binary Flow in Unfertilized Hamster Ova
Conclusions

Recomendations for Future Work

vi

vii

«hldld

11
15
25

33

42
50

56

63
63

66
69
73
78

81
81

93

103
111
118
118
131
137

141

Appendices:
Appendix A:
Appendix B:
Appendix C:
Appendix D:
Appendix E:
Appendix F:
Appendix G:
Appendix H:

Bibliography

Subprogram Unit MARBOX

Subprogram Unit MODEL

Subprogram Unit RK4

Subprogram Units DEQmD and CEXT
Subprogram Units TAKEPT and CIRCLE
Subprogram Unit CUHGEN

Tabulated Osmotic Shrinkage Data
Tabulated Binary Flow Data

iv

146
146
162
167
171
176
182
186
188

190

LIST OF TABLES

TABLE 5.1.1 - Results for Shabana's cell 1.
TABLE 5.1.2 - Results letting cod float as parameter.

TABLE 5.1.3 - Final permeability values as determined by parameter
estimation routine.

FIGURE 2.3.1
FIGURE 2.4.1
FIGURE 2 .5 .1
FIGURE 2.6.1
FIGURE 4 .1 .1
FIGURE 4 .1 .2
FIGURE 4 .1 .3
FIGURE 4 .2 .1
FIGURE 4.2.2

FIGURE 4.3 .
FIGURE 4.3 .

hihi

LIST OF FIGURES

The two-chamber system.

Membrane cross-section.

Membrane with capillary structure.
Shape of V. as a function of b.

Cross-section of the diffusion chamber.

Exploded view of the diffusion chamber.

Permeability of Cuprophan dialysis membrane as a
function of molecular weight.

Tho chamber system used to characterize permeability
of the dialysis membrane.

Typical results of dialysis membrane permeability
experiment for glycerol.

Typical ovum photographs from diffusion chamber.
Typical liposome photographs from diffusion chamber.

Point groups selection sequence and integration region
selection sequence.

cell
cell
cell
cell

response of
response of
response of
response of

and predicted
and predicted
and predicted
and predicted

Measured
Measured
Measured
Measured

‘WNH
I.

and 6.
and 8.

ova 5
ova 7

response of
response of

and predicted
and predicted

Measured
Measured

vi

NH" <4l480m-2$5~lf§°o OMU>

a min e.n-m'= «ruin

dc
a.

NOMENCLATURE LISTTNG

Area; membrane surface area.
Non-dimensional constant in Section 2.6.
Vector of parameter estimates in Chapter 3.
Local concentration (oamol/kg).
Concentration (oamol/kg).

Diffusion coefficient in free solution.
Proportionality constant.

Distribution coefficient.
Phenomenological coefficient.

Hydraulic or solvent permeability.
Number of moles.

Modified permeability coefficient.
Resistance phenomenological coefficient.
Sum of the squares function.

Time .

Temperature. absolute.

Velocity.

Specific volume.

Volume. .

Spatial coordinate.

Arbitrary driving potential.

Measured data values.

Greek Symbols
Non-dimensional constant in Section 2.6.
Vector of parameter values.
Model estimates to function values.
Chemical potential.
Osmotic pressure.
Volume fraction.
Reflection coefficient.
Variance of data point (i) in Chapter 3.
Entropy dissipation function.
Solute permeability.
Dialysis membrane permeability.
Non-dimensional time.

‘vii

CHAPTER 1

Introduction

1.1 Motivation £25 Prgsggt Iggk

The work described in this thesis was performed at the
Bio-Engineering Transport Processes (BTP) laboratory of Michigan State
University. The problem addressed was to adapt and develop the exper-
imental and analytical tools needed to perform experiments to
determine the passive transport preperties associated with the mem-
brane of an individual cell. This work is part of an ongoing study of
the effects of cryopreservetion procedures on the survival of dif-

ferent cell types.

The transport properties of a cell will have a direct effect on
the response of that cell to a freezing procedure. As water freezes
the solutes in the water are excluded from the crystal structure of
the ice. Thus as an ice front moves towards a cell the solutes are
concentrated into the remaining liquid water. This rise in solute
concentration will cause the cell to respond. in the form of an
induced volume change. in an attempt to regain an equilibrium state
with respect to the chemical potential inside and outside the cell.
The dynemic response for such a situation will be governed by the per-

meability characteristics of the cell.

 

During the freezing process two major mechanisms are thought
destructive. causing life-threatening damage to a cell. [22]. If the
rate of freezing is "slow" then the cell is subjected to high concen-
tration extracellular solutions for a relatively long period of time.
The cell will generally respond by expelling water from inside the
cell thereby decreasing the cell volume. In a sense the cell is
attempting to increase the intracellular concentration in order to
reestablish equilibrium. At slow freezing rates so much water is lost
that the cell can suffer damage due to the high intracellular and/or

extracellular concentration of solutes.

On the other hand. if the freezing rate is "fast" then the cell
will contain a large fraction of the initial amount of internal water
after extracellular ice has formed. This situation results in a high

prohability that internal ice ‘ill form. [23].

From this argument it is plain that the rate at which a cell
transports water and other permeable solutes will. to a great extent.
determine the freezing rate at which the cell is most likely to sur-
vive. It would be desirable to be able to predict the optimum
freezing rates before attempting a freezing process. This requires
that one know the transport properties of the cells of interest. It
would also require that the properties be determined as a function of

temperature as well.

The techniques described in this thesis make it possible to

determine the passive transport properties of an individual cell.

This is a significant advance over previously published methods (such
as the stop-flow methods) which only give average data based on bulk
samples of the cells. [2.8.12]. Being able to observe individual
cells as they undergo a non-equilibrium passive transport process
makes it possible to determine information on the distribution of
individual parameters within the population. It also makes it possi-
ble to work with cells for which a large bulk sample is not readily

available.

The diffusion chamber microscope stage described in Chapter 4.
Section 1 of the present work. which is used to gather the experimen-
tal data on the cells. can easily be modified to include a heat
exchanger system. This will allow one to study the effects of tem-
perature on the passive transport properties as well. This feature
was not yet installed at the time this thesis was completed so that
all of the data presented in the present work was gathered at room

temperature. .

1.2 Statement 21_Objgctizgs

The overall objective of this work was to adapt existing and
develop new experimental and analytical tools required for the deter-
mination of cell membrane permeability coefficients. It is intended
that this thesis will serve as a starting point for future investiga-

tors working on the problem of passive membrane transport.

The first step towards realization of this goal was to conduct a
literature survey to define the state of the art with respect to the
most popular and widely used of the passive transport models. The
presentation of these models is not intended to be all-encompassing or
fully detailed. Rather they are meant to be used as an introduction
to each of the models presented and to the field of passive membrane
transport in general. Once the reader has gained a basic understand-
ing of the passive transport process and the approach taken by various
investigators. it should be possible to work with the original presen-
tations of the models reviewed in this work and others not included

here with the fundamentals in hand.

The most widely used of the passive transport models. developed
by ledem and Iatchalsky in 1958. has been coded into FORTRAN for com-
puter implementation. This subprOgram is described in detail and
listed in Appendix C of the present work. This subprogram is intended
not only as a usable routine but also as an example of the interfacing
required between the modeling subprogram and the other subprogram

units described in this thesis. From this example the user should be

5
able to code and implement other models of interest for use with the

other subprograms developed as a part of the present work.

The actual determination of the values of the permeability param-
eters contained in the model of interest is accomplished through a
parameter estimation routine. Several of the available parameter
estimation routines are presented in chapter 3 of this thesis. Here
again the presentation is intended to serve only as an introduction to
the science of parameter estimation. By reading this chapter the user
should be able to gain sufficient understanding of the methods of
parameter estimation to be able to use the subprograms deve10ped as a
part of the present work. and to understand some of Ithe pitfalls
involved. In parameter estimation the pitfalls are many and the cor-
rective actions required are often learned only through experience.
The programs presented in this work have been extensively documented
and thoroughly tested and are to the best knowlege of the author in

good working order.

The parameter estimation subprogram MARBOX is presented in Appen-
dix A of this thesis. This subprogram is actually three parameter
estimation routines in one. It has the capability to run as the Ordi-
nary Least Squares (OLS) method (see chapter 3 sections 2 and 3). or
the Box-Kanemasu method (see chapter 3 section 4). or as Merquardt's
method (see chapter 3 section 5). It also enables the user to specify'
upper and lower bounds on the values of each of the floating parame-

ters.

 

The user is encouraged to make his own copies of this routine and
to modify it to suit his own needs. This will be particularily
appropriate in the adaptation of new models and in tailoring the input
and output formats. Each subroutine is written in a stuctured form in
order to facilitate understanding. Emphasis has been placed on docur
menting each routine with liberal use of comment statements and

variable definition blocks. The user is encouraged to follow suit.

The data reported in this thesis was collected using a simple
diffusion chamber microscope stage developed by Ligon and documented
in an unpublished work. [19]. This chamber is described in Chapter 4.
Section 1 of the present work (see figures 4.1.1 through 4.1.3). This
chamber makes it possible to subject an isolated sample of cells in
suspension to a psuedo-step change in extracellular concentration.
The cells remain stationary and are not sheared during this process.
The response of an individual cell to this treatment can be observed
directly. and videotaped or photographed. The resulting photo images
can be processed to yield cell volume as a function of time. This
processing can be done either by hand or through the use of computer
image analysis techniques. The computer is able to define and enhance

the boundary of the cell in the photcmicrograph.

This volume-time history. along with the experimental conditions
and initial estimates of the parameter values. form the input to the
parameter estimation routine. The parameter estimation routine then
calculates the statistically-determined values of the permeability

parameters which result in the best fit between the experimental data

 

7

and the values predicted by the model.

The overall routine has been tested using several sets of data
for two cases of the Kedem and Katohalsky permeability model. The
first case was that of osmotic shrinkage in which no permeable solute
is present. Osmotic shrinkage data reported by Shabana. [20]. for
unfertilized hamster ova was processed using the parameter estimation
routine linked to the Redem and Ratchalsky permeability model for this
case. The results are consistent with those generated by Shabana
using~ a- closed-form approximate solution to the Kedem and Katchalsky
equations. Data has also been generated on the response of unfertil-
ized hamster ova in a binary flow situation. that is. one in which
there is a permeable solute present. Processing of these data using
the parmmeter estimation routine also yielded consistent results.

The data on the response of the eve in the binary flow situation
was generated from photomicrographs taken with 35mm black and white
film. These photos were processed using simple computer image ana-
lysis techniques to yield the radius of the ovum as a function of

time.

Thken together these techniques form a powerful and versatile
tool for the determination of the passive transport preperties associ-

ated with the membrane of an individual cell.

CHAPTER 2

Historical Development of Permeability Models

2.1 Higgggical nggyigg

Characterizing the flow of materials across the membrane of a
cell has long been a problem of concern to investigators. Eerly
models. such as Jacob's model. 1952. [1]. attempted to describe the
passive transport process for membranes with expressions similar to
Fick's Law describing free diffusion. In this model the flowrate of a
given species was assumed to be linearly dependent on the spatial gra-
dient in concentration for that species. While this model adequately
described the transport process for some cases it was found that it

did not hold true in general.

The groundwork for the models generally used today was laid by
Onsager in 1931. [4]. Omsagor extended Lord Rayleigh's. [24]. work to
include thermodynamic flows and forces. Rayleigh originally expressed
a linear dependence between all mechanical flows in a system and all
the mechanical forces acting on the system. Thus Onsager proposed
what are referred to as the phenomenological equations to govern thor-~
modynamic systems which are not too far removed from equilibrium. He
also derived certain restrictions on the resulting coefficients based

on statistical considerations.

 

9
In 1958 Redem and Katohalsky. [3]. used an Onsagor set of pho-

nomenological equations and applied the principles of irreversible
thermodynamics in order to derive a model for the passive transport
process in a membrane. The resulting set of equations resolved the
inadequacy of previous models and has become the classical model for
membrane transport. This set of equations is generally referred to as

the REX permeability formulation.

The Kr! model was expanded.somowhat by Papanek in 1978. [12].
Papanek was able to derive a more general set of equations of the same
form as the EEK equations but without assuming ideal. dilute solu-
tions. The value of this aspect of Papanek's work remains to be
proven as the increase in generality comes at the expense of consider-
able computational effort and little comparative data between the
simpler I?! model and the Papanek.model is available. The Papanek
model will not be used in the present investigation as the considera-
tion of non-ideal solutions requires the determination of certain
empirical relationships for the solutes of interest and this is beyond
the scope of the present work. Papanek's model is presented in anti-

cipation of future work to be conducted in passive membrane transport.

Other investigators have derived direct closed-form functional
volume-time solutions for certain special cases of the KER equations.
[13.14]. Johnson and Wilson. [8]. derived an approximate closed-form-
power series solution for the binary KER equations. Direct solutions
of this type are generally easier to evaluate as they do not involve

the evaluation of a set of differential equations over time as in the

10
full set of Redom and Katohalsky equations. While these closed form

solutions may involve simplifiing assumptions and hence may not be as
accurate as the full set of Kedem and Katchalsky equations it may be
possible to use the simpler solutions to provide beginning estimates
to the values of the permeability coefficients at a relatively small
computational effort. This hypothesis has not been orplorod in the
present work but again the methods are presented in anticipation of

future work.

 

2.2 lggobg' gags;

One of the earliest models for the transport of materials across
a membrane was proposed by Jacobs in 1952. [1]. He considered the
problem of two regions in thermal and mechanical equilibrium separated
by a semi-permeable membrane with nonreloctrolyto solutions of differ-

ing concentrations on each side of the membrane.

The flow of solute through the membrane was described by an

equation analogous to Fick's Law:

dNildt = isA(c3-c§) (2.2.1)

where the superscript (0) represents one region and (i) the other; the
subscript (s) implies solute; N represents the number of moles; A the
membrane area; es the concentration of solute (in.moles/liter); and k3
is the proportionality constant or solute permeability coefficient

with units of (cm/sec).

Similarly the volume flow. which Jacobs related to water flow. is
assumed to be proportional to the transmembrane difference in chemical
potential of the solution. Assuming no hydrostatic pressure differ-
ence exists across the membrane. the chemical potential difference is
equivalent to the osmotic pressure difference. (xi-n")1 . Further
assuming ideal solution behaviour in both regions one can use the

relation’:

 

1 Iatchalsky and Curran. [9]. eq(10-8). pg.118 .
2 Iatchalsky and Curran. [9]. oq(5-55). pg.56 .

11

12

n = RTe (2.2.2)

where c is the osmotic concentration. which is the sum of the concen-
trations of all solutes. The osmotic concentration can be expressed
as the sum of the concentration of the permeable solute and the total

impermeable solute concentration:

0 = e, + c,ll (2.2.3)

The subscript (m) implies the sum of all impermeable solutes and the
subscript (s) implies the permeable solute. Thus the equation des-

cribing volume flow can be written:

dVi/dt = k,A(c1-c°) , (2.2.4)

where Vi is the internal cell volume. and k" is the permeability

coefficient for water ( the factor RT has been absorbed into k' giving

it typical units of cm‘lmole-sec ).

For a system of only one permeable solute one can write:

dVi/dt = k'a[[(ug+N§)/vi] - (cg+cg)] (2.2.5)

This is the form of the equation generally referred to as Jacobs'~

equation.

Many investigators have shown this model to be inadequate to des-

13
cribe the general passive transport process. Among those are Zeuthen
and Prescott. [2]. In their investigation Zeuthen and Prescott sub-
jected frog eggs (otherwise in equilibrium with the solutes present)
to high external concentrations of heavy water. 9,0. The heavy ygter
was shown to act like any other solute as its penetration into the
eggs followed the response predicted by (2.2.1) exactly. However.
while the solute penetrated the cell as expected. it was observed that
the volume of the cell remained constant so that (dV/dt = 0). From

(2.2.3) we see that:
dVildt = o = k'A(c1-c°) - (2.2.6)
' which implies that c1 = co. Using equation (2.2.4) then:
c: + c: a cg + cg (2.2.7)
As stated above the eggs were brought to equilibrium with all other
solutes present prior to immersion in the heavy water solution by

first immersing them in a similar solution with normal water replacing

the heavy water so:

°3 = c: (2.2.3)
This leads to the conclusion that:
cg - c: (2.2.9)

 

14
and this is clearly a contradiction of the known experimental condi-
tions. This then proved that (2.2.4) and (2.2.5) are incomplete

descriptions of the volume transport process.

Zeuthen and Prescott went on to show that for cells not in equi-
librium, with the non-penetrating solutes the penetration of heavy
water was not adequately described by (2.2.1). Penetration was found
to be more rapid in solutions where the extracellular concentration of
nonrpermeating solutes was lower than the internal concentrations and
slower in the reverse case. This demonstrates the need for some form
of coupling between the volume and solute flows which is not accounted

for in the Jacobs model.

2.3 The fizz Pgrmeability Formulation

In 1958 Kedom and Hatchalsky put forth what has come to be a
classic model for membrane permeation based on the principles of irre-
versible thermodynamics. [3]. This set of equations is known as the

RH! formulation and is still used extensively in various forms.

The EEK model resolved the inadequacy of the Jacobs model by
relating the £10108. Ji' of each species to all of the driving poten-

tials. X1. in the system through coupling coefficients. Lij‘

I1 3 LIIXI. + L131: + e e e + 111an
I: - L31X1 + L331: + e e e + 143an

. (2.3.1)

In 8 Ln1x, + Ln3X, + . . . + Lnan

These equations also can be expressed in their conjugate form where
the driving forces are expressed as functions of the species fluxes.
This leads to the so-called resistance formulation and will be dis-
cussed in the following section. It should be noted that use of these
equations implies the assumption that the system is not too far
removed from equilibrium as the linear relationship between forces and

fluxes can not necessarily be expected to hold as the forces increase.‘

The development begins by considering a system of two chambers

separated by a membrane. Each chamber contains a solution of a single

15

16

permeating solute in the same solvent (presumably water). The regions
are assumed to be in thermal equilibrium characterized by a single

temperature. T (see Figure 2.3.1).

 

r—Membrone

Chamber 0

 

 

 

 

 

 

Figure 2.3.1. The two-chamber system.

One of the basic principles of irreversible thermodynamics is
that if an adiabatic system undergoes a change of state via a reversi-
ble process then the entropy of the system will remain unchanged. If.
however. the adiabatic process is irreversible then the system will
experience a not increase in entropy. The rate at which entropy is
produced in a system which is in thermal equilibrium can be expressed
as the sum of the product of the flows in the system and their cor-
responding driving forces. The system of interest is defined as the

chamber designated chamber (i). which will be referred to as the cell.

and includes the membrane itself. no rate of entropy production for-

this system is given by the expression‘:

 

1 Kedem and Ketchalsky. [3]. eq(12).

17
dis/dt = (1/r) [(ue-u$)dN$/dt + (pg-p§)dN§/dt] (2.3.2)

where S is the entropy of the system. I is the absolute temperature. u
18 the chemical Potential. N} is the number of moles of component (j)
inside the cell. and t is time. The grouping di()/dt implies internal
to the system. The superscript (i) implies inside the cell while the
superscript (0) implies outside the cell. A dissipation function per

unit area of membrane. A. is defined for convenience as:
n - (TIA)diS/dt = (1/A)[(perui)dNé/dt+(u:-u§)dnildt] (2.3.3)
letting each mole flux be represented by: 55 = (1/A)dN}/dt (2.3.4)
o - ("S-vi". + (pg-phi, (2.3.5)

Thus for this lumped analysis the dissipation function. 0. is the sum
of the products of the fluxes h and their corresponding forces (the

differences in chemical potentials).

The system of particular interest in this study is that of a mem-
brane separating two solutions each made up of many solutes of which
the membrane is only permeable to one. Further. we are interested in
the case where there is no transmubrane hydrostatic pressure differ-
ence (AP-0) as most biological membranes will not support such a

difference. [3].

An alternate set of fluxes will be defined for use in the follow—

18

ing developent. The total volume flux through the membrane is given

by the expression:

The diffusion flux through the membrane is given as:

I. - (5,5,) - (5.1...) (2.3.1)

“10!. t1“ quantity 3, is defined by the relation:

AO‘IE' II 13(02/01‘) (2.3 .8)

If the transmembrane concentration difference is small then Ac‘/3'<(1
and 3,2(o§+c2)/2. Kodem and Ratchalsky describe the diffusion flow as
the "relative velocity of solute versus solvent which is a measure for
exchange flow". [3]. Note that c'-(1-§,-§1)/;' where i: represents the
volume fraction and the subscript (s) implies the permeable solute and
(i) the sum of all the impermeable solutes. Thus if the solution is
assumed dilute so that the total volume fraction of the solutes is
3.411 compared to 1 then 6,-1l—v'. This allows one to express Jd as:

Id = (is/3.) - in. (2.3.9)

In order to derive an expression for the entropy dissipation

function in terms of these new fluxes we must define a set of conju-

 

l9

gate forces for this set of flows. Let the conjugate forces for Jv
and Id be represented by TV and 1d respectively. Using these fluxes
and their conjugate forces the dissipation function analogous to

(2.3.5) then becomes:
0 ' (ifiaifiaxv + [(£,IE,)-(?,i,) 11d (2.3.10)

Because the entropy dissipation function.must remain unchanged under
the transform one can equate the two expressions for 0 represented by
(2.3.5) and (2.3.10) to form a single expression which excludes the
term 0. Th0 £13808 5. and n, are independent so their coefficients on
each side of this newly formed equation can be equated yielding two
independent expressions relating the new forces. X, and Id. to the old
set of forces. Au. and An'. These two expressions can be solved for
IV and In in terms of Au' and An'. The resulting expressions are

given by:

xv - o'au' + 3,1)”, (2.3.11)

x‘l - (l-fi.);'Au. - {‘O‘Ap‘ (2.3.12)

where {.e;;;,. If the solutions are assumed to be ideal and the

volume fractions of all the solutes are assumed to be small then’:

Au, - -(1/o,) [RTAc‘ + RTAci] (2.3.13)

2 Kedem and Katehalsky. [3]. eq. 35.

20
where the subscript (i) implies the sum of the impermeable solutes and

the subscript (s) implies the permeable solute. Similarily for the

permeable solute:

An. = arAc,fEs (2.3.14)

Introducing (2.3.13) and (2.3.14) into (2.3.11) and (2.3.12):

IV = -RTAci (2.3.15)

Id 3 RTACS + g‘RTACi (2.3 .16)

Using the Onsagor phenomenological equations each of the fluxes
in the system is assumed to be a function of all of the driving poton- .
tials in the system. For the case of two permeating species (the
permeable solute and the solvent) one will have two independent flows

and hence two independent forces related by the expressions:

II a L11x1 + Lizxs

The (L's) in these expressions are called the phenomenological
coefficients and are governed by Onsager's Law which requires that the

cross coefficients be equal. [4]:

L13 3 L31 (2.3.18)

21

For the system under consideration here (2.3 .12) becomes:

J'v ‘- I"va + I"pdxd

The values of the coefficients LP’ Lpd’ and La are restricted by the
requirement. under the principles of irreversible thermodynamics. that
the entropy production and hence the dissipation function must be
greater than or equal to zero. Substituting (2.3.19) into (2.3.10)
using (2.3.6) and (2.3.7) the dissipation function can be expressed

a a prv’ + “1.33de + ded’ 2 0 (2.3.20)

Since either Iv or Xd can be made to go to zero independently this
restricts both LP and Ld to positive values only. and requires that

the magnitude of Lpd be such that:

Most investigators will use a transform changing from the phe-
nomenological coefficients LP’ Ld. and I‘pd to an alternate set of

coefficients L1” 0‘. and a). [3.6.7.8.9.10]. The Staverman reflection

coefficient. a. is defined by the relation. [6]:

0' B [1 + All";';‘Aus]ngo (2 e3 .22)

 

22

The condition 0f vao occurs when the solute and solvent are flowing
in Opposite directions with magnitudes such that the volume of the
cell remains constant with time. Using (2.3.6). (2.3.7). and (2.3.19)

one can show that:

The solute permeability coefficient. 0. is defined by the relation.

[12]:

is 3 [ 633A“, 11 =0 (2.3.24)
v

so that using (2.3.6). (2.3.7). and (2.3.19) it can be shown that:

a) - ;‘[([.pLd - Hip/LP] - (Ld - Lpa’fi, (2.3.25)

The restrictions on LP and Ld and that represented by equation
(2.3.17) imply a restriction on m such that «20. The solvent permea-

bility coefficient. Lb' remains unchanged with respect to the

transformation.

By introducing (2.3.23) and (2.3.25) into (2.3.19) one can show

that:

Iv ‘ ’I-pRTEAn - deli“. + téAci] (2.3.26)

Id - mg. + «(Tu-an)“. + arm/3,) + «1.911., (2.3.21)

23

Since one must usually keep track of the internal solute content in
order to calculate the internal solute concentration it is often con-
venient to work with the total volume flux. Jg, and the solute mole
flux. fi,. rather than Jv and Id. From (2.3.6) and (2.3.7) one can

show that:
I'18 a (JV + Id);s (2.3 .28)
so that:
1is ' ;e(1-‘)Jv + uRTfAc‘ + Cgifici] (2.3.29)
Ignoring the contribution of the term {'EAci, which represents
the contribution of the impermeable solutes to the force la, in each
expression (as it is normally very small compared to Acs)' these

expressions reduce to the more commonly used set:

IV = -LPRT§Aci - oLpRTAcs (2.3.30)

.1, .. -u;,[Lp(1-.)§A.i + [oLpfl-o) - (u/E,))Ac,] (2.3.31)

It may be of interest to express the solvent flux alone. Using

(2.3.6) and (2.3.7) one can show that:

 

3 Kedem and Katchalsky. [3]. pg. 238.

24

i, - -(RT7;;)[Lr2Ac1 + (LPG + v‘)Acs] (2.3.32)

It is worth noting once again the rostictions under which this
set of equations remain valid. The system considered was that of two
compartments separated by a membrane and in thermal equilibrium.
Ihile there may be many solutes present on either side of the membrane
the membrane is assumed to be permeable to only one of those solutes
and the solvent. Further. the solutions are assumed to be ideal and
dilute due to the approximations for the chemical potential used.
Also it has been assumed that the driving forces in the system are
sufficiently small such that a linear relation exists between all the
driving forces and the resulting flow of each species. For the final
version of the equations presented. (2.3.30) and (2.3.31). it is also
assumed that no hydrostatic pressure difference exists across the mem-
brane. This assumption is not inherent in the overall development so
that if a pressure difference is believed to exist in the problem of

interest this can be accounted for in the I?! model. For situations
which comply reasonably well with these restrictions the I?! formula-
tion has been found to provide quite reasonable correlation with

experimental data. [8.10.11].

 

2.4 1;; EEK Rgsistance Formulation

In a work published in 1961 Kedem and Kachalsky presented an
alternate formulation of their membrane permeation model. [7]. This
formulation loads to a set of frictional coefficients governing the
process rather than the phenomenological coefficients as derived in
the previous section. [3]. Those frictional coefficients allow for a
more direct physical interpretation of the permeability parameters LP'
o. and a. While this formulation leads to more complicated expres-
sions than the previous formulation it has the advantage of physical
interpretability. It also represents an advantage in that the mechan-
ical coefficients encountered are largely concentration independent
whereas the phenomenological coefficients are generally concentration
dependent. [7]. While this model may not be practical for calculation
purposes. and will not be used for such in the preseht investigation.
as it introduces additional unknown factors such as the distribution
coefficient of solute in the membrane. it does provide some interest-

ing insights into the passive transport process.

In this formulation relations reciprocal to those represented by
(2.3.1) are used. In this form the driving forces. xi, in the gygtem
are assumed to be linear functions of each of the fluxes. J1, in the

system: .

25

 

26
X1 = R11J1 + 3121: + . . . "' Ran-Tn

x3 = R3111 + Rzsz + e e e + Rszn

. (2.4.1)

X1. =- Rnng + RmJ, + . . . + 11“an

Here the coefficients Rij are in essence frictional coefficients. The
numerical values of these coefficients are restricted by the require-
ment that the entropy production and hence the dissipation function
represented by (2.3.2) and (2.3.3) must be positive. This restricts
the straight coefficients. R11, to positive values and the magnitude

of the cross coefficients. Rij' to the condition:

Rh 1 gun” (2.4.2)

These coefficients must also satisfy Omsager's Law so that. [4]:

As before. a system of one permeable solute and one solvent (usu-
ally water) will be considered. Again we will assume that the system
is composed of two chambers separated by a membrane and in thermal
equilibrium. The solutions in both chambers are assumed to be
well-stirred so that no unstirred layers exist at the boundaries of
the membrane. The membrane thickness will be represented by the quan-
tity Ax so that the membrane covers the region 0<x(Ax (see Figure

2.4.1).

27

The local molar flux of solute and solvent can be written as:

Is = Csvs = dnsldt

I' a: C'V' = dn'ldt (2.4.4)

where n1 is the number of moles of species (i) per unit area and Vi is

the local velocity of species (i).

r—\

Membrane
Solution 0 Solution i

.____.x

/"—\(

L—AX—-l

Figure 2.4.1. Membrane cross-section.

 

 

In order to deve10p a set of local equations which can be
integrated across the membrane the system will be assumed to be in a
steady state condition. This will allow us to express various quanti-
ties as total (derivatives with respect to x rather than having to
concern ourselves with a set of partial differential equations in x
and t. This implies that the concentrations. C1, are only functions
of x. as are the velocities. vi. Note also that the concentration.

C1. is defined as the number of moles of species (1) per unit volume
of membrane not per unit volume of solution. This dofinition.makes C1

a local concentration accounting for the membrane volume as well as

28
the volume of the solution. The fluxes. 1i: will be assumed to be

independent of x and equal to the macroscopic or total transmembrane

fluxes.

The problem is subject to the following boundary conditions:

us(x=0) = u:
u'(x=0) = u:
n,(x!Ax) - a:
H(x=Ax) - .13, (2.4.5)

Since the objective is to derive a local set of equations to be
integrated across the membrane. local driving potentials must be used.

The local driving force is the local gradient in chemical potential:

Xj 3 'dllj/dx (2.4.6)

(which has units of dyne-cm/mole-sec). This is consistent with the
development which led to equation (2.3.5) except that in section 2.3
the internal mechanics of the membrane were not considered and a set
of transmembrane equations. as opposed to a set of local equations.

was being developed. Using (2.4.6) for this system (2.4.1) becomes:

-d"w’d‘ ' 8waw + RwsIs
-du,/dx - 33,1, + 1:”1, (2.4.7)

where: Rs' 3 2"

29
It is now assumed that under steady flow conditions the thermody-

namic force acting on a given component will be counterbalanced by a
sum of mechanical friction forces so that there will be no net accel-
eration. These forces are assumed to be comprised of a set of forces
due to a given species interacting with each of the other permeating
species in the system and a force due to each species interacting with
the membrane itself. For the two component system being considered
here this translates to:

X. = -F - F

8' sm

x' a ’Fws — me (2.4.8)

where F represents the frictional force and the subscripts s. w. and m
refer to solute. solvent. and membrane respectively. Thus the factor

F‘h would represent the frictional force between the solute and mem-
brane per unit mole of solute. Each of the frictional forces is in
turn.represented by an expression relating the force. Fij' to the
relative velocity of component (i) to component (j) through a mechani-

cal friction coefficient. fij' per mole of (i) component:
Choosing the reference frame such that vat: 0. (2.4.6) becomes:

"dllsldx = f"(v,-v') + fmv8 = v8(f"+fm) - v'fs'

-dp‘ldx - f"(vi-vs) + f'nv' = “stws + v‘(f',+f'n) (2.4.10)

 

30
It is also possible to rewrite this set of equations as:

-du,/ax . [(rflumwcshs - (fs'IC'N,

-du,/ax = -(£,,/C,)J, + [(£,,+fm)/c,]1, (2.4.11)

Comparing (2.4.11) with (2.4.7) readily shows the frictional

coefficients to be:

3“ 3 (f3'+f‘.)/Cs 3" = -f,'/C'

a" .. -f'3/Cs a" = (f's-I-meC' (2.4.12)

The restriction represented by Onsager's Law requires that:

-fsw/cw = 'fws/Cs (2.4.13)

Thus only three independent frictional coefficients exist while the
fourth is dependent. This is analogous to the need for three indepen-

dent permeability coefficients in the model presented in section 2.3 .
[Odell and ntchllsky point out that "f" is of the “no mtnre as

the friction coefficient of free diffusion f:' given by Einstein's

equation:

if“, .. um' (2.4.14)

and may therefore be assumed to be approximately independent of the

31

local concentration as D° is essentially concentration independent".
Do is the diffusion coefficient of the solute in free solution. This

implies that due to the relation expressed by (2.4.13) the coefficient

f" must be strongly dependent on the solute concentration. This
would imply that substituting fsw for fws would be the most reasonable
strategy in order to yield friction coefficients which are independent

of the concentration. This leads to the expressions:

-du,/dx = [(f..+fm>/C.]J. - (rm/can,

’dl-l'ldx = -(£,,,/c‘,);rs + [[(cslgnflnmJ/CJJ, (2.4.15)
or:

'Csdusldx = (£8'+£sm)13 - (Cs/C')fs'J'

4'4“”, - ‘fst + [(cs/c'nnumh, (2.4.16)

It must be noted that (2.4.16) represents a set of local equa-
tions only and that in order to be truly useful as experimental tools
they must first be integrated across the membrane. The integration
procedure is presented in detail in Kedem and Katchalsky's paper. [7].
and for the sake of brevity only the results are presented here.

After integrating from 0 to Ax (2.4.16) becomes:

‘Ans ' ' w[fsw;wi;;s/§wm]Ax + J.(f"+f‘.)Ax

§n(AP-Aui-Ans) =- J,[£n+(£"?,I;F,)/§“]Ax - Isrwu (2.4.17)

32

where 5'“ is the volume fraction of water in the membrane. and the
quantity Kc°s represents the mean value of solute concentration across

the membrane and is given by:
.!_. Ax 1
co, = (l/Ax) L chsdx (2.4.18)

where K: is a local distribution coefficient of the solute in the mem-

brane given by:
I: - eye: (2.4.19)

Cs represents the local concentration of solute based upon the total
volume of the membrane itself. The quantity cs represent. the concen-
tration of a free solution of equal solute chemical potential.‘ If the
solution is ideal in both chambers then!c is constant and:

to 8 tides]. (2.4.20)

1:6, = ridufic‘; + c§)/2 (2.4.21)

A further simplification also occurs if the solute and solvent
penetrate the cell only by passing through solution-filled capillaries
in the membrane rather than through reaction with the membrane itself.
Iedem and latchalsky observed that for a case such as this the solué
tion in the capillaries approaches the composition of the free
solution. This means that the solute distribution coefficient becomes

equal to the volume fraction of water in the membrane.

2.5 Comparison 2; KFK Resistance in; Permeability Formulationg

By comparing the results of the two different formulations of
Kedem and Katehalsky presented in the previous two sections interest-
ing relationships appear. These relationships give added meaning to
the physical significance of the permeability parameters which are
commonly used to describe transport of materials across a membrane.
The relationships derived here will not be used for calculation pur-
poses in the present work.but rather are presented for the insights
into the passive transport process which they provide. An alternate
method of comparison to that presented by Kedem and Katchalsky. [7].
will be presented here. The overall strategy will be to pose the
resistance equations represented by (2.4.17) into the same form as the
permeability equations with flows expressed as functions of the driv-

ing forces.

Let us begin by recalling the EEK permeability equations in the

form:

I, - -(1/v')[LpAni + (LP. + 3,4011%] (2.5.1)

I8 = -';,[[Lp(1-a) - 39.1.1“ + [chU-c) - (./E,)]An,]_ (2.5.2)

where we have included the CsEhci term which was neglected in (2.3.30)
and (2.3.31). The comparable set of resistance equations. represented
by (2.4.17). now must be expressed in this same form. One way to do
this is to solve both of the equations for I, and equate the results.
The resulting expression can then be solved for J' to give:

33

34

Kfs‘Ans - §m(fs'+fsm)(A“1+Afls) (2.5.3)

 

Afom( fs'4—fm) + hfs'fm]

where: k -‘V'ch,/§'n. A is used only for convienience of presenta-
tion. Note that it has been assumed that no transmembrane difference
in hydrostatic pressure exists (AP=O) in order to be compatible with

the Kedem and Iatchalsky set.

In a similar manner. by solving both of the equations represented

by (2-4-17) for 1'. equating the results. and solving for Is it can be

shown that:

J’ ___ “fwnflfsflms ' tmlfMAnfiAn.) (2.5.4)
3

 

Axtf“(f"+fm) + “at,“

It is now possible to equate the coefficients of the independent
osmotic pressure differences to obtain the desired cross relation-
ships. For instance by equating the coefficients of Aui in (2.5.1)

and (2.5.3) it can be shown that:

2,2,. (2,, + 53,) (2.5.5)

 

Axlfn(f"+fa) + 128,1...)

A check on the units for this equation shows the right-hand side to
have dimensions of (length’lforce-time) which are the proper dimen-

sions for LP' the solvent permeability.

35

This equation also can be written:

1 2.5 .6

 

Lp = (E'gm/Ax)[fm + 1f."

In this form it is apparent that LD 1; invorgoly proportional to the
sum of the solvent-membrane friction factor. f'm, and an additional
factor characterizing solute-solvent interactions (the bracketed
term). The direct dependence of solvent permeability. Lb’ on volume
fraction of solvent in the membrane is also apparent. This seems con-
sistent as a membrane with a larger solution content would be expected
to be more permeable than an otherwise identical membrane with a

smaller solution content as the larger solution content would imply

that a larger volume fraction of the membrane is available for flow.

Kedem and Katchalsky developed an expression for Lp only for g
single. very limited case. It would be expected that (2.5.6) should
reduce to the I?! expression under the same restrictions. Tb show

that this is the case we begin by rearranging (2.5.6) as:

 

1 )TI- (2.5.7)

L =- " Ax f +2.: (
p (v.6...) )[wm sm “(fa/f")

Consider now only the case of a coarse non-selective membrane. Kedem

and Katehalsky give the non-selectivity condition as:

fem/:8 ' fn’;w (2.5.8)

36

using this in (2.5.7) one can show that:

 

- 7 - — 1 1 (2.5.9)
1.1) = (v.4.n/Ax)[fm[1 + Mvslv'H 1+(fu/f") )1]—

Kedem and Katchalsky further restrict the consideration to a case with
no transmembrane concentration difference. This reduces the value of
the mean concentration of solute within the membrane to the concentra-
tion of the free solution.multiplied by the (constant) distribution

coefficient:

2'5. " ‘5. (2.5 .10)

While there would be no flow in this case one can still derive a sim-
plified expression for the solvent permeability. Just because there
is no flow this does not imply that the solvent permeability does not

exist.

For membranes with a capillary structure (see Figure 2.5.1) the
distribution coefficient is equal to the volume fraction of water in

the membrane so that:

l = V'c' (2.5 .11)

and (2.5.9) becomes:

37

 

. .12
1 )Id (25 )

L = '- -'-
p (v'tm/Ax)[fm[l + vsc.( 1+(fsm/fsw)

Membrane

 

Capillary

Flow

Figure 2.5.1. lbmbrane with capillary structure.
Unless the membrane is somehow actively inducing the flow of either
the solvent or the solute the frictional coefficients will be posi-
tive. That is to say that the presence of the membrane will act to
inhibit the flow. For relatively dilute solutions with positive fric-

tion coefficients:

$33, << 1 (2.5.13)

1 + (fa/t") > 1 (2.5.14)
so that we can further reduce (2.5.12) to:

LP = (3.1:...) Han“) (2 .5 .15)

38
This is indeed the same expression derived by Kedem and Katchalsky for
Lp under these conditionsI. This lends greater confidence to the val-
idity of the more general expression for L represented by (2.5.6)

P
derived for the first time in this work.

Further confidence in the validity of (2.5.6) in particular and
to the present strategy as a whole is gained by consideration of the
Staverman reflection coefficient. By equating the coefficients of Ans

in (2.5.1) and (2.5.3) one can show that:

_. ._ If .-
a = -(v,u/Lp) - (Vw/LPA1)[ " g“'(fs'i'fsn) ] (2.5.16)

f'n(f"+fm)+lf“f8m

Substituing for LP using (2.5.5) in the last term only and rearranging

gfves: -
a =- 1 - (V‘u/Lp) - [(n") / [snuflnmn] (2.5.17)

This is again identical to the relationship derived by ledem and

Iatchalsky. This same result can be achieved by equating the coeffi-

cients of Ari in (2.5.2) and (2.5.4).

Kedmm and latchalsky point out that if the solute and solvent
penetrate the cell by different paths so that there is no

solute-solvent interaction then f"p0 . Using this in (2.5.17) gives

 

1 Kedem and Iatchalsky. [7]. eq. 4-21.

 

39

the expression:

a = 1 - (V‘s/LP) (2.5.18)

This relationship then can be used as a test for the condition of
non-interaction between solute and solvent flows. It is interesting
to note that while this appears to imply that the solute and solvent
flows are uncoupled they are in fact coupled in the sense that the
magnitude Of An, effects the solute flow as well as the solvent flow
and that the magnitude of Aus effects both the solvent flow and the
solute flow. This coupling is implied in the phenomenological equa-

tions.

In order to carry the deve10pment of o farther more information
about the nature of the solute permeability m'is required. One way to
solve for m as a function of the friction coefficients would be to
equate the remaining coefficients (those of An8 in (2.5.2) and
(2.5.4)). substitute for L1) and c. and solve for m. This however
requires solving an extrmnely complex quadratic equation. For this
reason only the form developed by Kedem and Katchalsky will be consi-

dered here.

Kedem and Katchalsky developed the following expression for the

solute permeability:

a - n1 - :,(1-.)1 I mu" + run (2.5.19)

40
In.most cases the solute volume fraction. :3, will b. voxy small com-
pared to 1.0. This would reduce the above expression to the more

simple form:

a) = x / [Ax(f" + {an (2.5.20)

If the membrane is assumed to be of the capillary structure then K=§wm

and this expression reduces further to:

w = 2:... I [Ax(f" + me (2.5.21)

These relationships show that the solute permeability is inversely
prOportional to the sum of the solute-water and solute-membrane fric-
tion factors.

It is possible to use the above relationships in (2.5.17) to
develop an expression for the reflection coefficient. or to use
(2.5.17) in (2.5.19) to develop an expression for solute permeability
in which only the friction coefficients appear. These procedures lead
to complicated expressions which are not particularly useful and hence

will not be presented.

In summary this comparison has shown that L1, can be oxprgsggd go
a function of all the friction factors with fwm plmyin‘ the donlngnt
role (see equations 2.5.12 and 2.5.15). The reflection coefficient.
a. was shown to be a function of the ratio of the solute permeability

to the solvent permeability. "/Lb' and an additional factor character-

.(I

41
izing solute-solvent interactions. An expression for (n was not
derived using this method but Kedem and Katchalsky showed that a) is
dependent on the factors f" and fun which characterize the

solute-solvent and solute-membrane interactions respectivly.

2.6 A nger Series Solution to; he §:§,Pgrmeabilitz Equations

In 1967 Johnson and Wilson. [8]. developed an approximate solu-
tion to the Kedem and Katchalsky permeability equations represented by
equations (2.3.30) and (2.3.31). The solution they presented was
based on a perturbation analysis and a power series expansion of the

K9! equations.

While the values obtained from this solution may not be as accu-
rate as one might wish they can be used as starting values for other
more accurate parameter estimation routines. The model presented here
gives volume as a function of time directly from a closed-form analyt-
ical expression and hence is readily evaluated. Other methods which
deal with the full set of K9! equations in differential form will
require the repeated numerical integration of the governing equations
and hence will require much greater computational effort. In these
latter types of routines good starting estimates of the parameters can
significantly reduce the number of iterations required to reach the
final values. In many cases this may mean the difference between a
routine converging to a solution or not. Thus by using the model of
Johnson and Wilson as the generator of starting values for other more
accurate routines one should be able to realize a significant reduc-
tion in total computational effort required to estimate permeability
parameters. While this method has not been utilized in the present

work it is presented in anticipation of future work to be conducted.

The development begins by considering the I?! equations in the

42

43

form:

awe: - LPRT‘Aljucon-tchHV/l - cos - c.] (2.6.1)

st/dt = earA[e, - (NS/W] + E,(1-.)dV/dt (2.6.2)

where N: is the number of moles of the permeating solute inside the
cell. V is the volume of the cell. V. is the initial cell volume. A is
the surface area of the cell. c, is the initial concentration of all

is the external concentration of the

the impermeable solutes. cs

P05301519 8013t°o and 3; is the mean tranmnmnbrane concentration of
the permeable solute defined by equation (2.3.8). LP’ a. and o are
the solvent permeability. solute permeability. and reflection coeffi-
cient respectively. These equations carry the implicit assumptions
that the cell is in equilibrium prior to time zero so that the
impermeable solute concentration is the same inside and out. and that
none of the permeating solute is present inside or outside the cell.
At time equal to zero the extracellular concentration of the permeat-
ing solute undergoes a step change from zero to c‘ and remains

constant thereafter. It will also be assumed that the surface area of

the cell remains constant.

Johnson and Wilson define an alternate set of permeability coef-

ficients by combining the factor RT with LI) and m such that:

p, - Lpn-r (2.6.3)
P - «RT (2.6.4)

44
This gives P; typical units of (cm‘lmole-sec) and P typical units of

(cm/sec). Using (2.6.3) and (2.6.4) in (2.6.1) and (2.6.2) one gets:

dV/dt = AP'[(o,v. + ensw" - arcs - c.] (2.6.5)

st/dt = pA[e, - (us/w] + gamma: (2.6.6)

This set of equations is then nondimensionalized using the following

groups:

v‘ = Vlv. N‘ = N8/c3V.
t . P'A¢.t/V. b = P/P'c.
u = cos/oo (2.6.9)

V7I is the volume nondimensionalized with respect to the initial
volume. N. is a nondimensional permeable solute content. and t is non-
dimensional time. The factors b and a are used for convenience of

notation. Using these groups in (2.6.5) and (2.6.6) yields the

expressions:

V'dV‘Idt = l - V. + c(N‘-V‘) (2.6.8)

v‘dN‘lde =- b(v‘-N‘) + (1-e)(3,/e,)v‘dv‘/dc (2.6.9)
This set of equations is subject to the initial conditions:

v‘(0) .. 1 N‘(0) = o

45

In a typical case the solvent will penetrate the membrane faster
than the solute so that the cell will initially shrink in size as the
solvent will be leaving the cell faster than the solute is entering.
Eventually the internal solution will reach a concentration high
enough to cause the solvent to begin reentering the cell and the cell
volume will increase. This is the typical "shrink-swell" behaviour
observed for many cases. At the point where the cell reaches its min-

imum volume dV‘ldr-O so that by equation (2.6.8) we see that:

V; = (1+uN;)/(l+u) (2.6.10)

where the subscript m implies minimum. Since N;)o one can show that:

V; > (1 - a) ‘ (2.6.11)

Thus for small values of c (which from (2.6.9) implies small changes
in the permeable solute concentration and/or little rejection of the
permeable solute) the perturbations in volume will be small. One can

then express V.(r) and N.(t) as power series in a such that:

v‘(.) - V:(:) + uV:(t) + a‘v:(c) + . . . (2.6.12)

N‘(c) - N:(t) + uN:(t) + e‘N:(e) + . . . (2.6.13)

where V;(t) and N;(:) are independent solutions which when weighted as
indicated by powers of u and summed up will yield the full solutions
V. and N‘. Note that for values of a less than one. un90 as n96.

Thus for small values of a a small number of the independent solutions

46
will have a significant contribution to the total solution.

The initial conditions on v'(t) and N‘(t) can be satisfied by the
power series expressions by requiring that v:(o)=1; v;(o)=o for (1)0);

and that N;(0)=0 for all (i).

These power series expressions can be substituted into (2.6.8)
and (2.6.9). One can then collect like powers of a in the resulting
expressions and equate the coefficients of a given power of u on
either side of a given equation. Equating the coefficients of the
zeroth power of a in the equation resulting from (2.6.8) gives the

expression:

V: dV:/dt = 1 - v: (2.6.14)

This equation and the initial condition on v: gra gatigfigd by v:=1.

Equating the coefficients of the zeroth power of c in the equa-

tion resulting from (2.6.9) gives the expression:
V: dN:/dt a b(V: - N:) + (1-a)(3./c,)v:(dv:ldt) (2.6.15)
Using V: - l in this expression reduces it to:

4N:ldt - b<1 - NI) (2.6.16)

This expression and the initial condition on N: are .otlgflod by:

47
N: = l - exp(-br) (2.6.17)

Equating the coefficients of the first power of a from (2.6.8) gives:
.

V:(dV:/dt) + v:(dv:/dt) = -vI + N: — v.

or: dV:/dt = v: - exp(-bt) (2.6.18)

This equation and the initial condition v:=o are gatisfigd by:

VI = -(1-b)-1 [em-ht) - mum] (2.6.19)

Ignoring the higher order terms one can now write:

v’ = v: + av: - 1 - [a/(l-b)][oxp(-bt) - 0*P“"] ‘2"°2°)

This expression will have an error of the order a3 which is

approximately of order (1-V:)’.
Johnson and 'ilson also give a similar expression for the case of
cells initially in equilibrium with a penetrating solute present and

subjected to a step change in extracellular solute concentration from

the initial value to zero at time equal to zero. This expression is:
V. -“-' 1 + [u'/(1-b)][exp(-bt) - exp(-1:)] (2.6.21)

where c'toN;(0)/c.V. and N,(0) is the total number of moles of perme-

48

able solute in the cell at time equal to zero. This case corresponds

to the ”washing” of a solute from a cell.

In order to obtain the desired permeability parameters from the
above model we need three pieces of information relating the experi-
mental data to the model. The first of these is the matching of the
shape of the curve for the experimental data to the shape of the curve
for the model. Jahnson and Wilson suggest using the following expres-

sion to characterize the shape of the modeling curve:

[exp(-bt) - exp(-1:)] (2.6.22)

. t
(V'-V n)/(1-V I!) = 1 " [Olp(’b“.) _ .xP(—1n)]

Note that this expression does not involve the variable a. Figure
(2.6.1) shows the shape of the curve described by (2.6.22) for.various
values of b. This figure shows that reciprocal values of b will give
the same curve. Johnson and Wilson observe that this implies that two
different sets of values for the permeability parameters will satisfy
the model equally well. They found however that one of these two sets
would include unrealistic values such as c>l and could therefore be
eliminated. They also observe that the model will be insensitive to b
for values of b=1 so that special care must be taken in a computer

program implementing this solution. '

Once the shape of the experimental curve has been.matched to that
of the model and a value for b has been obtained. one can then match

the time and volume scales of the model to that of the experimental

49
data. The time scale is matched by matching the time at which the

minimum volume is reached. From (2.6.20) one can show that:

t; = (ln b)/(b-1) (2.6.23)

80 that a value of t; can be calculated. Similarly one can match the
volume scale by matching the values of the minimum volume for the
model to that of the experimental data:

v; e 1 - ablb/(l'b)] (2.6.24)

and a value of a can be obtained. Using the numerical values of a. b.
and r and the definitions of u. b. and t given earlier one can deter-
mine the values of the permeability parameters LP’ w. and o as one

will have three equations and three unknowns.

 

 

 

 

'° I I I I I
an
00. 0-1}
3 06 0-0.5 1
F 6-20. 5'5
7 04» A
8
0-2— '“
l I I I I I
o l 2 5 4 5 e r a

Figure 2.6.1. Shape of V. as a function of b.

2.7 In; nggnek 319le

In 1978 Papanek. [12]. presented a permeability model which
closely paralleled the Kedem and Katchalsky permeability model. [3].
The Papanek.model differs from the Kedem and Katehalsky model in that
the assumption of dilute and ideal solutions is not made. Thus the
Papanek model attains broader applicability at the expense of computa-
tional simplicity. This model will not be used for calculation
purposes in the present work as it requires the deve10pment of certain
empirical relationships for each of the solutes of interest in order
to handle the non-ideality. This is beyond the scope of the present
work. A Since the develOpment parallels that of Kedem and
Katchalsky many references will be made to the equations of section 3

of chapter 2 in the present work.

Papanek begins by assuming that the phenomenological equations.
(2.3.1). hold true. The same expressions as those used by Redem and
Iatchalsky are also used for the local rate of entropy production.
(2.3.2). and the entrOpy dissipation function. (2.3.3). Papanek's
model differs from the Kedem and Katchalsky model in that Papanek does
not substitute for the difference in chemical potential. Au. but rath-
er retains it as the driving force. Kedem and Katchalsky on the other
hand assumed that the solutions are dilute and ideal and transformed
from the Au driving force yielding the hydrostatic pressure difference
and the difference in solute concentrations as the new driving forces.
By retaining the chemical potential as the driving potential Papanek

removes these restrictions from the model.

50

51
Papanek defines a simple set of phenomenological equations for

the situation involving a binary flow of water. w. and a single per-
meating solute. s. using the molar mass flux. defined by (2.3.4). and
the difference in chemical potential (Ap=u°-u1) as the flows and

forces respectively:

a O I
n' = J" 3 inAl‘w + L”Al-ls

'5. = I. = LLAIL. + L:2Alls (2.7.1)

where Onsager's law. [4]. requires that L:,=L:1. Substituting these
relationships into the dissipation function (2.3.5) one can show that:

0 = 1,1“, + ISA“, (2.7.2)

0 = 1:1,“: + 21;,AHM. + LL”; (2.7.3)

The fact that the rate of entrapy production and hence the dissipation
function must be positive definite implies the following restrictions

on the phenomenological coefficients:

.
L11 2 o
. ,
L33 1 o (2e7e4)

<LI.)‘ 1 LLLL

.Next the alternate set of flows derived by ledem and Iatchalsky are

umed:

52

v - 75,, + 13?. (2.7.5)

2...
l

where Iv is the total volume flux and Id is the velocity of solute
relative to solvent. The phenomenological equations for this set of

flows are:

Iv = 14.x. + Lidia (2.7.7)

Id L;1Xv + L;,xd (2.7.8)

where Onsager's law requires that L;,=L;1. The corresponding dissipa-

tion function for this set of flows and forces is:

It is now required that'the local rate of entropy production and hence
the dissipation function.must remain unchanged by the transform. This
implies that one can equate the two expressions for 0 ((2.7.2) and

(2.7.9)). Recognizing that 3' and Is are independent one can then
equate their coefficients using (2.7.2). (2.7.5). and (2.7.9) to show

that:

X. = Au'l(;'(1+0)) + Lats/(1+0) (2.7.10)

xd = [A6, - (?./?,)A..'] 3,!(1+e) (2.7.11)

where 06?;33. Using these expressions in (2.7.7). (2.7.8). and

53
(2.7.6) one can show that:

c...
ll

,, [$.(1w)’]"[L;.l(AH/3.05.2163] + L;,[;,Aus(l-O)-29Au'l-vw]

' LLIOQAII. - anal-v.1] (2.7.12)

. ['c'./(1+e)’l[L;,[(A..,/?,)+E,Apsl + L;.[((1-o)hu,/?,)+2E,Aus]

H
ll

+ 1%. [F.An.-(ehu,/?,) 1] (2 .7 .13)

A conversion can now be made from the phenomenological coeffi-
cients to the more widely used set of permeability coefficients LP’ a.
and c. The first of these. LD’ is equal to the first of the phe-
nomenological coefficients. L;1. This parameter is normally called
the solvent (or water) permeability and has typical units of

cm'ldyne-sec. The solute permeability at zero volume flow. m. is

defined by the relationship:

J8 - [Jun] J =0 (2.7.14)
v
The situation of Jv=0 occurs when the solute and solvent are flowing

in opposite directions and with magnitudes such that the volume of the

cell remains constant with time so that:

43,5, - 3,1,, (2.7 .15)

54
By eliminating Au, from (2.7.13) using (2.7.10) and (2.7.5) one can

show that:
m = r8./(1+e)1[L;.L;.-(L;.)’I / [L;.-0L;.I (2.7.16)
Note that a has typical units of mol/dyne-sec.

The final parameter. c. is the solute reflection coefficient

which is defined by the expression:
a= 1+1 ” 2347
[ u'/(v'c,Au,)]Jv§0 ( )
so that using (2.7.13). (2.7.10). and (2.7.5) one can show that:

a = -L;.(1+e) I (L;,-OL;.) (2.7.15)

Substituting these new parameters into (2.7.12) and (2.7.13)

gives:

1' - [$..]"[[LP/e +«e;,.16e,/?; + [Lp(l-o)/a - $.613.Au,] (2.7.19)

I. a .“[IL,(1-e)/. -'3.eIE.Ae,/?;.+ [3.Lp(1-e)'/a + 913.4u.] (2.7.20)

where s-l+O-o0. Papanek points out that the factor a is typically

very close to one and hence can usually be considered to have very

55
little effect on the evaluation of the above expressions. Papanek

prOposes as a worst case an ethylene glycol solution with 38=4M.
I;,=42.3cm’/mol. 6-0.8 so that for this highly concentrated solution

s=1.03=1.0 so that one can. in general. assume s=1.0.

While this set of permeability equations is less restrictive than
the Kbdem and Katchalsky set one must deal with the problem of calcur
lsting the chemical potential of both the permeating solute and
solvent inside and outside the cell. This adds to the computational
complexity of the model. Since very little comparative data between
the simpler I?! model and the Phpanek model is available it is still
unclear whether or not a comparable increase in accuracy is also

achieved.

2.8 Noggimensiggglization of thg EEK Peggeabilitz Eggations

For practical calculations it is desirable to work with a nondi-
mensional (or partly nondimensional) set of equations. In this
section the Kedem and Katehalsky equations for coupled binary flow
represented by equations (2.3.30) and (2.3.31) will be nondimensional-
ized. The equations used for calculation purposes in the parameter
estimation routines deve10ped as a part of the present work are the
partially nondimensional set represented by equations (2.8.14) and
(2.8.15). The fully nondimensional set of equations are not used as
one must absorb either the solute or solvent permeability into the
nondimensional time factor to achieve full nondimensionalization.
Since the permeabilities are unknowns in the parameter routine this
makes the independent variable time an unknown factor as well. This
unneccessarily complicates the input of data to the routine. Thus the
set of equations used for practical calculations retain the dimension-
al time factor as the independent variable. Only the cell volume and

cell solute content are nondimensional.

Begin by writing (2.3.30) and (2.3.31) in the form:

JV . "and“iwi’ + a(c:-c§)] (2.8.1)

6, - shrug-ob + amen. (2.3.2)

where the superscript (0) implies outside the cell and (1) inside the
cell. Reacall that these equations apply to cases where there is no

hydrostatic pressure difference. AP-O. and where the solute volume
56

57
fraction is small. The subscript (i) implies the sum of all imperme-

able solutes and (s) the permeable solute.
Recall that Iv is the volume flux and can be expressed as:
IV = (1/A) dV/dt (2.8.3)
and that h. is the permeable solute flux given by:
h. a (1/A) dN/dt (2.8.4)

where A is the surface area of the membrane. V is the cell volume. N
is the number of moles of permeable solute inside the cell. and t is
time. An alternate set of permeability coefficients will be used and

are given by:

P, a LPRT'
(2.8.5)

Pa = uRT

so that P, '111 have typical units of (cm‘lsec-mole) and P. .111 have

typical units of (cm/sec).
Using (2.8.3). (2.8.4). and (2.8.5) in (2.8.1) and (2.8.2) gives:

avm - -P,1[(eg-ei) + e(c:-c’;)] (2.8.6)

dN/dt - P.A(O:-c:) + Z.(1—e)dwdt (2.8.7)

58
If it is assumed that there is a relatively small difference in the

concentration of the permeable solute across the membrane then the
mean permeable solute concentration. 33. which is defined by equation

(2.3.8). can be approximated by the arithmetic mean:
3, a (camp/2 (2.8.8)

Also notice that the internal concentrations can be expressed as func-
tions of the cell volume. The internal permeable solute concentration
is equal to the number of moles of permeable solute inside the cell.

N. divided by the cell volume. V:

c: = N/V (2.8.9)

Similarly the impermeable solute concentration inside the cell can be

expressed as:
“i " c1.V./V (2.8.10)

where °io is the initial impermeable solute concentration and V. is

the initial cell volume.
_Incorporating this information into (2.8.6) and (2.8.7) gives:

dV/dt - -P,I[(.g - chm/v) + .(c‘; - NIV)] (2.8.11)

ledt - P.A(c: - N/V) + (112m): + N/V)(1-c)dV/dt (2.8.12)

59

Defining the nondimensional volume and solute content as:

v‘ = V/V. N’ = NVc°v. (2.8.13)

and using them in (2.8.11) and (2.8.12) gives:

dv‘ldt = —(P,A/v.)[(eg - ei./v’) + ecg(1 - N‘Iv‘)] (2.8.14)

dN‘Idt - (P,A/v.)(1 - n‘lv‘) + (1/2)(1 + N‘Iv‘)(1—e)dv‘/dt (2.8.15)

While equations (2.8.14) and (2.8.15) are not fully nondimensional-
ized. as they still contain the time factor. they are a good set of
equations to use for practical calculation. This is because the data
available is typically the nondimensional volume as a function of real
time a priori, In order to fully nondimensionalize the equation set
one must absorb either P, or B. into the time factor. Since P8 and P'
are both typically unknowns this makes it impossible to calculate the
nondimensional time as a function of the real time. This will make
the procedure for estimating the permeability coefficients unnecessar-

ily complicated.

If one _wishes to work with a fully nondimensional set of

equations then the groupings used by Johnson and Wilson. [8]:

t a P'Aci.tlv.
b - P'IP'ci. (2.8.16)

a . ace/c1.

60
can be used so that (2.8.14) and (2.8.15) become:

dv‘ldt = [(1/v‘) - (cg/81.)] + a[(N‘/v‘) - 1] (2.8.17)

dN‘Idt = b[1 ~ (N‘lv‘)] + (1/2)[1 + (N‘Iv‘)](1-e) dv‘ldt (2.8.18)

For the case of no impermeable solute present at the initial state

(°i930) the groupings in (2.8.16) will not work. For this case one

can use the alternate set of nondimensional groups:

1' = PQAcgt/V.
b' = P,/P,cg (2.8.19)

a' - «we:
in which case the analogous nondimensional equations become:

av‘xdtv . -[.v - (ch/.gv‘) + an - (N‘Iv‘)1] (2.8.20)

dN’Idt' - b'[1 - (N‘Iv‘)] + (1/2)[1 + (N‘Iv’)](1-e) dv‘ldt' (2.8.21)

No further complication occurs if no impermeable solute is introduced
at time equal to zero. One can merely set cgso in (2.8.20) and

(2.8.21).

Mhny authors will simplify the evaluation of the above equations
by assuming that the volume changes induced in a cell will be suffi-
ciently small such that the surface area of the membrane. A. remains

constant [10.12]. The surface area is then typically calculated using.

61
the initial cell diameter. Stusnick. [10.11]. explored the effects of

this assumption on the resulting values calculated for the permeablil-
ity coefficients using data reported by Stewart and Jacobs. [18].
Stusnick first calculated the permeability coefficients assuming that
the surface area remained constant and then.by allowing the surface
area to be calculated as a function of cell volume using the expres-

sion:
A = (36n)"” v’“ (2.8.22)

While for the most part the change in the results was not significant

the parameter values in at least one case varied by as much as 20%.

Papanek points out that for a typical case of a cell subjected to
a rise in the extrhcellular concentration of.the permeable solute
while in the presense of impermeable solutes a shrink-swell behavior
is observed: "The rate of rapid initial shrinkage will be determined
almost wholly by the water permeability. LP' a. LP))m and there 1;
initially very little solute transport. As solute begins to enter.
and water continues to exit the cell. the interaction term a becomes
important. The shape of the volume minimum seen ........ will be
fixed primarily by c. As solute slowly enters. the rate of swelling
is almost entirely a function of solute permeability a. During this
phase. intracellular concentration is increasing and water re-enters
the cell. but the change in cell volume is rate-limited by permeable

solute entry”. [12].

62

One would expect these trends to express themselves in reverse in
a parameter estimation routine. For instance the value of LD will be
determined almost wholly by the rate of initial shrinkage. In the
case of the shrink-swell behavior the assumption of surface area
rmuaining constant at its initial value will consistently over-predict
the actual surface area. This would suggest that permeability would
be correspondingly under-predicted in order to keep the volume flux
the same for an over-predicted area. particularly during the initial
period of rapid shrinkage. This effect does appear when comparing the
values of LP in Stusnick's work. The effect does not. however.
express itself in the values of w whose value was changed in alternate
directions under the new evalustion depending on the specific case.
It is also interesting to note that the value of c was consistently
over-predicted when using the assunption of constant area as compared
to the values when the area was calculated as a function of volume.
One would suspect that this might be due to compensation for the

under-predicted values of Lb' .

With the availability of fast. efficient computers the extra
evaluations required for calculating the surface area of the membrane
as a function of volume will not significantly increase the execution
time and hence should. in general. be included in order to prevent

biasing of the calculated parameter values.

CHAPTER 3

Parameter Estimation lethods

LIMAi—tieal! ma 912117.191

When one thinks about solving an equation they inherently think
in terms of solving for the state of a system given certain parameters
and initial and/or boundary conditions. The parameters are given no
special consideration but are merely looked up in a hand book. in
many cases. however. the parameters are unknown or perhaps unknown
functions of the state variables. It is this problem which is

addressed by the methods of parameter estimation.

Pbrameter estimation is the science of determining from measured
data the values of the physical properties of materials and systems or
arbitrarily chosen parameters which play the role of constants in an
equation or mathematical model. The determination of physical prepar-
ties is typically performed through measurements of the state
variables (temperature. pressure. etc.). These measurements are
matched to a mathematical model of the process being monitored _in
which the parameters of interest appear either explicitly or implicit-
ly. In the case of arbitrary constants one may be trying to fit a
physical prOperty of a system to an empirical or analytical model.
This model may be a function of one or more of the state variables.

63

64
Typically one would be working from measured values of the preperty of

interest at specific states of the system. Many other applications
for parameter estimation methods have been proposed including the

separate yet related science called function estimation by Beck. [15].

It is not the purpose here to enumerate the applications of
parameter estimation but rather to present for the student new to
estimation theory a simple explanation of some of the more common
methods. This presentation is particularily slanted towards the prob-
lem of parameter estimation in membrane transport. The methods
applicable to this problem are not the most powerful techniques avail-
able. Because of the nature of the data being processed very little
information of a statistical nature is available. For instance. since
the data available using the techniques described in this work are
typically the volume-time history of an individual cell in a single
experiment no knowlege of the variance-covariance pr0perties of the
data is available. Given this type of information many more powerful
techniques could be employed. One of the simpler of the methods which
utilize statistical information. maximum likelihood (IL) estimation.
will be presented in order to illustrate one of the ways in which the

investigators knowlege can be incorporated into these methods.

The first methods to be presented are the ordinary least squares
(0L8) method and the Gauss method of minimization for nonlinear esti-
mation. These are perhaps the oldest and most basic of the parameter
estimation methods. Tbgether they form the basis of a multitude of

other techniques. Some of these modified least squares methods will

65
also be presented including the Box-Ianemasu method and larquardt's

method.

A FORTRAN computer program which is based on the OLS estimator
and the Gauss minimization method has been coded as a part of the
present work. This program also incorporates the modifications of
Iarquardt and Box-Kanemasu. This program is call lARBOX and is des-
cribed in detail in Appendix A of the present work. This program has
been particularly tailored to the passive membrane transport problem

which has at most three parameters.

66

3.2 Ordinary Least §gpp£gg pad lgxlmum ngglihgog Egtlggtegs

The Ordinary Least Squares (0L8) method is perhaps the oldest and
simplest of the parameter estimation routines. This method. like all
other estimation methods. attempts to minimize the error between the
function values predicted by a mathematical model and the known or
measured functional values with respect to the floating parameters.
The various routines will differ in the way in which they measure the
fit between the data and the model and how they choose the direction
and size of the step-changes taken to determine the "true" parameter
values from initial estimates. It is this method which has been
incorporated into the parameter estimation routines developed as a

part of the present work.

In the OLS method the fit between the data and the model is meas-
ured by the sum of the square of the difference between the measured
values and the corresponding values predicted by the model. This is

referred to as the sum of the squares function and for OLS is given

by:

S - Y (11 - MUD)a (3.2.1)
i-l
where n is the number of data points. the 11'. are the measured data
values. and the “1'8 are the values predicted by the model. ‘3
represents the vector of parameter values so that the notation n1(3)
expresses the dependency of the predicted values on the parameter

values.

67

The problem then becomes that of minimizing the functional 8 with
respect to the parameter values. lathematically this is expressed
through the derivatives of the functional S with respect to the vector

of parameters by meeting the condition:
asmr‘ - 0' (3.2.2)

where 5. is the vector of "best fit” parameter values. When the
values of F' are found that satisfy equation (3.2.2) then the sum of
the squares function has been minimized with respect to the parame-
ters. Since one can not generally satisfy (3.2.2) exactly one
typically will continue iteration on the parameter values until the
change in parameter values becomes insignificant. It is possible in
some cases to have more than one set of values satisfy this condition
due to the presence of local minimums. One would wish to find the
values corresponding to the global minimum. This generally is not a

problem as most well-posed models will possess only one minimum value.

Under laximum Likelihood (IL) estimation a slightly modified sum
of the squares function is utilized. Each data point is weighted by

its variance so that 8 becomes:

8 - 1 - ("))/ 1' (3.2.3)
In. 2—1 H 1 “19 ‘1

where G, here represents the variance of data point (i). Thus these
data points with a lower variance will be weighted more than those

with higher variances. This type of method is not utilized in the

68

present work as the varience of the data points used in the present
work are unknown. Mhny other methods which utilize the investigator's

knowlege are available.

3 .3 33 EM Minimization Mpihod

The Gauss or Gauss-Newton method for minimizing the sum of the
squares function is one of the simplest minimization methods. While
it was prOposed many decades ago it is only with the advent of the
computer that it has become practical to apply due to the large number
of calculations required. This method has proven to be effective for

well-defined problems which have a distinct minimum.

Recall that we are attempting to find values of the vector F such
that an equation of the form (3.2.2) is satisfied. By differentiating

(3.2.1) one can show that for the OLS method:

11
asla - a a - “ ‘ . .
3 I3 [2131 (11 “(an ] (3 3 1)
as/a a -2 Y - ‘)) a a" ' (3.3.2)
B 2:31 ( i 1|1(B “i/ B

Because 3 is a vector. satisfying (3.2.2) implies simultaneously

satisfying the set of p equations:

aSIap: a 0
68/88: a 0

. (3.3.3)

38/33; . 0

where p is the number of floating parameters.

69

70
A sensitivity matrix X is defined with the components:

xij = ani/apj (3.3.4)

so that X has dimensions of nxp. Further using the matrix notation on

Y and 11:

T
Y = [Y1.Y,......Yn] (3.3.5)
_ T
n = [al.n.......nn] (3.3.6)
equation (3.3.2) can be rewritten:
T __
88/01! = -2 x (Y - 11(3)) (3.3.8)
Thus the desired solution. 3.. will satisfy the condition:
T — a
X (Y-n('6 )) =0 (3.3.9)

The problem is now to solve for the value of 3‘ in this expression.

This is not easy,however,as both X and a are dependent on the value of

F.

In the Gauss method the first two terms in a Thylor series expan-
sion of i about the estimated values of 3.. designated 5. are used to
generate a better estimate of 3.. This method requires that (K have
continuous first derivatives in ‘3 and bounded higher order derivi-

tives. The Taylor series is given by:

71
fia‘) = ah) + (mm/ah) (ii’- 8') + H.O.T. (3.3.10)

or using the sensitivity matrix and ignoring the higher order terms

(H.0.T.):
- '6) (3.3.11)
Substituting this expression into (3.3.9):
T _._ _w. ._o ._ -
X [Y - nu.) - X(b)(B - b)] a: 0 (3.3.12)

This expression can now be solved for the approximate value of B. giv-

ing:
- r - 1' ..
II‘ :8 + IX X]1 x (Y- “(6n (3.3.13)

For models which are linear in the parameters. that is for models in
which the parameters appear only to the first power. equations
(3.3.11) through (3.3.13) will be true equalities and only a single
iteration will be required to reach the final values. However for
models nonlinear in the parameters. several iterations may be required
depending on the nature of the model and the accuracy of the initial

parameter estimates.

The iteration procedure can be summarized as follows. Defining

the matrix P as:

72
T -1.
P = [X 'X] (3.3.14)

and using the superscript (k) to indicate the current values and (k+1)
to indicate the future values then one iterates on the following two

equations:

1 T k k
P( ) = [‘2 ( )‘X( ) 1‘1 (3.3.15)

T _
8m” = 6”“ + 15ml 2 (“a - “mu (3.3.16)

where fi‘k) is a function of 5(k). The iteration continues until there
is a negligible change in the parameter values. One way of testing
for convergence is to set a tolerence factor. 8. such that execution
ends when:

ABSI bit“) _ his) ]

(k) < 6 (i=1.2...p) (3.3.17)
ABS[ bi ] + 81

 

where 51 is a very small number (81<<<8) which will prevent a divide
by zero conditon in the event that one or more of the parameter values
goes to zero during an iteration. Typically the value of 61 will be

about 10"16 but this will depend on the computer being used.

3.4 flpgzgggggggp Interpolatiop Mpihod

The Box-Kanemasu.method is a minimization method based upon a
modification of the Gauss method. In some cases this method can pre-
vent oscillations of increasing amplitude in the parameter values
caused by the linearization of a nonlinear model using the Gauss
method. Included in this presentation is a modification to the
Box-Kanemasu method used by Bard. [17]. which insures a decrease in

the sum of the squares function for each iteration.

The basic strategy in this method is to modify the step size cal-
culated by the Gauss method in order to insure a better fit with the
new parameter values. The direction of the step is not changed under

this method.

Begin by decomposing equation (3.3.16) into two equations:

A5 (7- 11

k (k _T k) _(k)
‘ ’= p ’[x( )] (3.4.1)

_ 1 _.k
b(k+ ) = F.(k) + h(k+1)Ab( )

(3.4.2)

where h is a scalar quantity whose magnitude can vary with each itera-
tion. If h-l for all values of (k) then one has the Gauss method.
Equation (3.4.2) is typical of many of the modified least squares
methods. These methods will differ in the way in which the value of h
is calculated. Sane methods will perform a search on h to determine
the value which gives the minimum value for the sum of the squares

function.
73

74
In the Box-Kanemasu method the value of h which will minimize the

sum of the squares function is approximated through a second order

quadratic approximation for S in h:

S = a. + alh + alh3 (3.4.3)

where a.. a1. and a, are arbitrary constants to be calculated for each
iteration. A corresponding linear approximation of the parameter vec-

tor 5 will also be used:

F=b+hfl (man

To calculate the values of the three arbitrary constants in (3.4.3)
three pieces of information are required. The first is to calculate S

at h-O so that:

(k)

S(h-0) = a, = s. (3.4.5)
_ k
'3(hp0) = b( ) (3.4.6)
1 _(k)
e. = s: ) a S(b ) (3.4.7)

A second value is determined at has where a is initially set to 1.0.
This value of S is denoted Six). The presence of u in this develop-
ment represents the modification used by Bard. i If c-l then the
following development would match that of Box-Kanemasu. In the modi-
fied version the value of e is manipulated to insure a reduction in

the sum of the squares function.

75

The final piece of information is gained from the derivative of S

(k)
with respect to h in the Ab direction at h=0:

01’:

From equation (3.3.8) for

[08/01.]th0 = e. = §:=1

- T
a1 = [(aSIaa) (dF/dh)]h=o

the

ordinary least

squares function one can show that:

T _
08/08 a -2 X (Y - n)

and from (3.4.4) one can show that:

so that:

(
The scalar quantity G

so that:

(k

dB/dh a 85

“=[4iTW-fiflTmu)

G(k) a [

_T

X

‘1 a —ZG

)

k
) is defined as:

- k
a-anm‘°

(k)

[33/331]h,0 [asi/ah]h=0

squares

(3.4.8)

(3.4.9)

of the

(3.4.10)

(3.4.11)

(3.4.12)

(3 e4 e13)

(3.4.14)

76
Using the definition of the matrix F given by equation (3.3.14) it is

possible to show that:

- k T ‘- _ k
a‘k) = [Ab( )1 ‘P 1 Ab( ) (3.4.15)
It can also be shown that:
k k k
a2 = [5; ) - s: ) + 26‘ )ulluz (3.4.16)

The minimum of S with respect to b can be found by setting the
derivative of S with respect to h equal to zero and solving the
resulting expression for h. This will yield the expression:

(k+1) (k) (k) (k) (k)
h = + 26

a u’/[Sa - so a] (3.4.17)

Bard suggests using the variable a to insure that:

s(k) (k)

a . (3.4.18)

If this condition is not met by the initial value of c=l then u is
51Vidfid by t'0 and 3a is recalculated. This process continues until a
value of u is found such that the condition of (3.4.18) is met. If
the value of o.must be reduced to about 0.01 and (3.4.18) is still not
satisfied then it is suggested that execution be discontinued. Beck.
[16]. suggests that this condition may indicate improper programming.

especially in the calculation of the sensitivity coefficients. or near

77

linear dependence of the sensitivity coefficients.

The use of equation (3.4.17) in the calculation of h is restrict-
ed by the value of 1,. Beck. [16]. explores the implications of all
three possibilities on the value of a,, that is. a,<0. a,=0. and a.)0.
The results show that equation (3.4.17) should be used only if a,)o.
This can be tested by requiring that:

(k)

k k
S( ) E ) - (2-A—1)uG (3.4.19)

0. is

where A is a value equal to or slightly greater than one. Beck sug-
gests using Ael.l. If the inequality represented by (3.4.19) is not
satisfied while (3.4.18) is satisfied then Beck suggests that:

k 1 '
h( + ) = An. (3.4.20)

k+1
be used. If (3.4.19) is satisfied then h( ) should be set to the

value calculated by (3.4.17) or by (3.4.20) whichever is the smaller.

3-5.!sraasrd£;s.!eihed

Marquardt's method of minimization utilizes a compromise between
two different methods. One method is called the method of steepest
descent. In this method the step direction is taken in the direction
opposite the gradient. Thus the direction chosen for the parameter
step is that of the steepest downward slope in the sum of the squares
function with respect to the parameter values. Mhthematically this is

expressed by:
Ab “
j ‘ -aslabj = 22131 (Y1 - ni)ani/abj (3.5.1)

This method typically works well when the parameter estimates are far
from the minimum but is slow to converge near the minimum as the value

of BS/abj approaches zero and the step size decreases.

The second method is the Gauss method of minimization which util-
izes a linearization of the fitting function to calculate the step
size. This method is outlined in section 3.2 of this work. The Gauss
method works well for parameter estimates close to the minimum but

poorly when the estimates are far removed and the model is non-linear.

' In Isrquardt's method a compromise is made between the two
methods. Initially the method resembles that of the steepest descent
and as the parameter estimates improve it shifts toward the Gauss
method. This strategy is implimented by weighting of the Iii matrix.

which is the matrix of sensitivity coefficients.

78

79
The matrix A is defined as:

A = X ‘X (3.5.2)

so that I is a pxp symetric matrix. Also note that from equation

(3.3.14):
P=x (3.503)
The matrix A is modified using the expressions:

A“ = Anne.) (i=j)
(3.5.4)

so that if i=0 the matrix is unchanged. For large values of 1 this
procedure weights the A matrix towards its main diagonal driving the
routine towards the steepest descent method. As k90 the method resem-

bles the Gauss method.

In application the initial value of l is made relatively large
(say about 1000). The worse the initial parameter estimates and the
more ill-behaved the model the larger the initial value of h. A new
set of parameter values and the sum of the squares function for these
new values are calculated based on the initial parameter estimates.
If the sum of the squares function has been reduced by the new parame-

ter values then A is reduced by a certain factor (anywhere from 2 to

80

about 10) and iteration proceeds. If the sum of the squares function
is not reduced then.l is increased by the same factor and a new set of
parameter values is calculated using this new value of A. This con-
tinues until a reduction in the sum of the squares is realized. The
routine then proceeds to the next iteration. New parameter values are
calculated based on the updated estimates and iteration continues

until a standard convergence criterion is met.

CHAPTER 4

Data Gathering Techniques

4.1 I}; Difiugion Chamber

The diffusion chamber is a microscope stage capable of subjecting
an isolated sample of cells to a psuedo-step change in extra-cellular
concentration of solutes. These solutes can be either permeable or
impermeable to the cell membrane. Figure 4.1.1 shows a cross section
of the chamber and Figure 4.1.2 shows an exploded view of the chamber

parts.

The sample of cells to be tested is held in a region bounded by
glass on top. a dialysis membrane on the bottom. and a rigid support
ring around the side. The thickness of this region is manipulated by
a spacing ring which is thick enough so that the cells are not com-
pressed between the membrane and the glass and yet not so large as to
allow the cells to drift in and out of focus. Thus the size of the
cells of interest determines the thickness of the spacing ring and

hence the sample region.

Below the sample region on the other side of the dialysis mem-
brane is a second region containing a bulk flow of the extracellular
solution. This region has a volume much greater than that of the sam-

81

”\O

MICROSCOPE
OBJECTIVE

 

 

 

GLASS\
RIGID DIALYSIS
SUPPORTT‘ SAMPLE 4 / MEMBRANE

J I l

DIFFUSION

 

 

 

 

 

 

 

 

 

 

 

 

—> BULK FLOW —>-

 

 

'1 fl 1
L I

LGLIIss

\ /
LIGHT SOURCE

Figure 4.1.1- Cross-section of the diffusion chamber.

82

View Hole - ‘—T0p Flttlng
U
@"— Spacer Ring

Dialysis
Membrane —.

 

 

 

 

 

Membrane
Retainer
Bulk Flow
Outlet Port
Bulk Flow ; .
Inlet Port 2
Plastic -
Bulk Flow Brass

 

 

 

 

Channel
\ 77

Glass Base

/ /

Figure 4.1.2- Exploded view of the diffusion chamber.

 

 

83

84
ple region. Initially the bulk flow region will contain the same

solution as the sample region. The chamber and the cells are allowed
to come to equilibrium with this solution. At time equal to zero the
bulk flow region is flushed with a new solution whose solute concen-
tration differs from that of the initial solution. If the new
concentration is higher than that of the initial solution then solutes
will diffuse upward across the dialysis membrane from the bulk flow
region to the sample region. If the new solution has a lower solute

concentration then the solute flow will be reversed.

The bulk flow can be regulated through a pair of pressurized bot-
tles. one containing the initial solution and the other containing the
new solution. The bottles are connected with plastic tubing to a pair
of small pumps with a single outlet. his outlet is connected to the
bulk flow inlet port on the diffusion chamber. A plastic tube is con-
nected to the bulk flow outlet port and inserted into a discharge
vessel. One switch is provided to switch the pump intake from one
bottle to the other and a second switch turns the pumps on and off.
This mechanism was developed by Ligon and is detailed in an unpub-
lished work. This technique works well for very large cells which do

not have a tendency to move about.

For smaller cells the pressurized bottles created large scale
disturbance due to the sudden introduction of flow in the bulk region~
which could cause the cells to drift out of view. For cases where
this is a problem an alternate injection method has been developed.

Rather than using the pressurized bottles a large (35cc) syringe is

85

used to introduce the bulk flow. The chamber is set up initially with
the equilibrium solution in the bulk flow channel. The syringe is
filled with the new solution and a plastic tube attached to the end.
Air is then removed from both the syringe body and the tube. The
other end of the tube is then connected to the bulk flow inlet port of
the diffusion chamber. At time equal to zero the syringe is manually
depressed introducing a bulk flow of the new solution into the bulk
flow region. This system allows the user to maniputate the velocity
of the bulk flow in order to prevent large scale disturbances in the
system. When using a syringe with smooth plunger action this techni-

que is easily implemented.

Once the fluid in the bulk flow region has been replaced by the
new solution the change in solute concentration will be gradually
introduced into the sample region. The cells in the sample region
will begin to respond to the non-equilibrium osmotic conditions creat-
ed. The response of the cells will depend explicitly on the
extracellular concentration of the sample region as expressed in equa-
tions (2.3.30) and (2.3.31). Thus in order to estimate the value of
the 90:3035111t7 P‘fletOI' Lb' c. and a one must have explicit
knowlege Of the °:(t) and c2(t). the extracellular concentrations of
permeable and impermeable (relative to the cells membrane) solutes.
These extracellular concentrations correspond to the concentrations in
the sample region so that the extracellular concentration.will depend
on the rate at which the new solution in the bulk flow region mani-
fests itself in the sample region. This will in turn depend on the

permeability of the dialysis membrane to each of the solutes present.

 

86

In Ligon's work it was shown that for most solutes the sample
region will have very small concentration gradients and hence can be
characterized by a single concentration. In the following a model is
developed which allows for the prediction of cs(t) and °i(t) in the
sample region based on the solute permeability characteristics of the

dialysis membrane used in the diffusion chamber.

The diffusion chamber itself can be modeled using the Kedem and
Ketchalsky definition of solute permeability. The defining relation-

ship for solute permeability. u. is given by equation (2.3.24):

53 . [ne,An,]Jv=o (4.1 .1)

Iwhere 3’ is the mole flux. m is the solute permeability. and 3' is the
mean' transmembrane solute concentration defined by equation (2.3.8).
A“, is the transmembrane difference in chemical potential in this case
given by ("Sb-“80) where (b) implies the bulk solution and (e) the
extracellular or sample solution. It will be assumed that the volume
of the sample region remains constant and Jv=0° In oxporimgnts run at
high magnifiction using the chamber the dialysis membrane does not
appear to deform to accommodate volume changes as the cells will gen-
erally remain in focus throughout the experiment. This would tend to
support the validity of the constant volume assumption. This means

that for each solute present one can write that:

Ji a 61318.11 (4.1 .2)

87

Using the expression. [9]: Afll - :3AP + RTAci/zi (4.1.3)

and assuming no transmembrane pressure difference (AP-O) then:

Again since the dialysis membrane does not appear to deform during an
experiment. and since the fluids are incompressable the assumption of
no hydrostatic pressure difference appears justified as no other
mechanism exsists by which a pressure difference could be maintained.
Defining N as the number of moles of a given solute in the sample

region and dropping the subscript (i) one can write:
dN/dt - JA.- eRTAAc (4.1.5)

where A is the effective transfer area determined by the spacing ring.
One can then use Ae'cb-c. and c‘-N/V where V is the constant volume of
the sample region. Note that if the cells make up a significant por-
tion of the sample region then V can be replaced by (V-Vc) where V0 is
the total volume of cells. If the overall change in cell volume is
assumed to be small then V can be assumed to be constant and no furth-
er change in the following development is needed. Equation (4.1.5)

then becomes:

dN/dt - mRTA [ob - (N/VI] (4.1.6).

or: dN/dt + (uRTA/V)N’- eRTAcb (4.1.7)

Assuming that °b is not a function of time (i.e. that it undergoes a

88

step change at time equal to zero) then this is a linear first order
ordinary differential equation in N and t. It is subject to the ini-

tial condition:

N(0) = N0 = cov (4.1.8)

This differential equation can be solved by two linearly
independent solutions which together satisfy the initial condition and
the differential equation. One solution is generated by considering

the homogeneous equation:

th/dt + («IRTA/V)Nh = 0 (4.1.9)

0

This equation is satisfied by the expression:

Nh = B exp[-(uRTA/V)t] (4.1.10)

where B is an arbitrary constant. The second solution comes from the

non-houogeneous equation (4.1.7) and is given by:
Nn a °bv (4.1.11)
Reconstructing the full solution:

N a thn a B exp[-(mRTA/V) t] + Vcb (4.1.12)

89
Applying the initial condition:

N. = B + Vcb (4.1.13)

.6 that: B = N. - v6b (4.1.14)

This makes the final solution:

N(t) = Vcb + (No-vcb)exp[-(mRThlV)t] (4.1.15)

Dividing this expression through by V and letting c(t)=N(t)/V then:

c(t) = 6b + (c, 7 ch).xp[-(naIaIV)t] (4.1.16)

Thus the concentration of each solute in the sample region is
expressed as an exponential function of the concentration of the bulk
solution. the initial concentration of the sample region. the permea-
bility of the dialysis membrane to that solute. the volume of the

sample region. the effective area of transfer. and time.

The most difficult of these factors to obtain will be the permea-
bility of the dialysis membrane. Data has been obtained for the
permeability of Cuprophan dialysis membrane produced by Enka
Glanzstoff AG. These data imply a linear relationship between the log
of the permeability of the dialysis membrane for a given species and
the log of the molecular weight of that species. This dependency is

shown in Figure (4.1.3) for three types of the Cuprophan.membrane.

91

It is interesting to note that the time constant for equation

(4.1.16) is given by the expression:

t. - V/AmRT (4.1 .17)

Consider as one example the typical configuration used to generate the
data reported in the present work. In that configuration the quantity
VIA was given as 2.286x10'zcm. The Cuprophan 80pm membrane was used
so that for sodium chloride «RThl.166x10"cm/sec. This yields a time
constant of 19.6 sec. Thus it will take 58.8 seconds (three time con-
stants) to reach 955 responce in the sample region when using sodium
chloride. This is a considerable length of time for experiments which

typically last 4-6 minutes.

Consider as a second example the same configuration of the diffu-
sion chamber and the same dialysis membrane but with sucrose as the
solute of interest. The molecular weight of sucrose is 342.3 grams so
that the permeability of the 80pm dialysis membrane to sucrose is
approximatly 3.3x10" and the time constant becomes tc-69.3 see. This

means that it will take nearly 4 minutes to achieve 95% responce.

From these simple calculations it is clear that the presence of
the dialysis membrane in the diffusion chamber will introduce a signi-
ficant delay in the time in which the sample region comes to
equilibrium with the new solution in the bulk flow region. One would
expect this delay to have a significant effect on the values of the

permeability parameters calculated in a parameter estimation routine

92
based on data gathered using the diffusion chamber system.

4.2 Characterization g; the Dialzgig Mbmbrane Peggeabiliiy

A technique deveIOped by Ligon. [19]. and modified as a part of
the present work has been used to characterize the permeability of the
dialysis membrane to various solutes. In the present study the Cupro-
phan 100pm flat membrane was tested using the solute glycerol. Other
studies performed in the BTP lab by Dupuis have utilized other solutes
including sodium chloride and sucrose. [25]. The values for the per-
meability of the dialysis membrane generated as a part of the present
work were used in the processing of data on the binary flow response
of hamster embryos. The results of this process are documented in
Chapter 5. Section 2 of the present work. The method used to charac-
terize the dialysis membrane permeability is presented here in order
to document the procedure used for future reference. The characteri-
zation of the dialysis membrane permeability to various solutes will
be an ongoing research project in the BTP lab. In processing data
obtained using the diffusion chamber system the dialysis membrane per-
meability value used in the transport model will have a great effect
on the resulting cell membrane permeabilities calculated as shown in
Chapter 5 of the present work. Thus it is important that the dialysis
membrane permeability to each of the solutes of interest 'be known
accurately in order to obtain reliable results using the diffusion

chamber system.

This technique utilizes two well-stirred chambers separated by a
piece of the dialysis membrane (see Figure 4.2.1). One chamber con-

tains 400ml of a relatively low concentration solution (typically

93

 

/Router
§ H

Ring stand

 

 

 

 

/

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

( 1
Stir rodsm‘~“| J/yb/
F
r ‘lzInner chamber
' Outer chamber
t__, w
Membrane
. {Retainer ./Stir-bar
Ring stand
m/z”/’ ii:j
” ‘ \ N
.1?— 73h Magnetic stirrer

 

 

 

 

E— .

Figure 4.2.1- Experimental set-up for determination
of dialysis membrane permeability.

94

95
0.1-0.2 moles/liter) of glycerol. The other chamber holds 1600m1 of

distilled water.

The larger outer chamber used in the present work was a rectangu-
lar (9x4x6 inches) chamber made from clear plexiglass and sealed with
silicon caulk. This chamber was set off-center on a magnetic stirrer
platform so that the stir-bar was positioned towards one end of the
chamber. The smaller inner chamber. which holds the higher concentra-
tion solution. was made from a section of 3 inch inside diameter PVC
drainage pipe approximately 8 inches long. A clear plastic window was
cut into the side of this chamber in order to facilitate observation
of the inner solution during an experiment. The lower end of the
inner chamber was milled to accommodate a membrane retaining ring also

made from plastic PVC pipe.

A piece of the dialysis membrane approximately 4 inches square
was soaked in distilled water for 20-30 minutes prior to use in an
experiment. At the beginning of an experiment the membrane was
stretched over the lower end of the inner chamber and the membrane
retaining ring pressed over the membrane. holding it securely in
place. The inner chamber was then suspended within the outer chamber
approximately 1 inch above the bottan of the outer chamber using a
ringstand and clamp. The actual level at which the inner chamber was
suspended was predetermined such that when the inner chamber held
400ml and the outer chamber held 1600ml the fluid levels in each
chamber would be the same. This insured that no hydrostatic pressure

difference would be created due to unequal levels in the chambers.

96

The presence of the window in the inner chamber made it possible to

check the levels visually as well.

The inner chamber was positioned at the end of the outer chamber
opposite the end with the magnetic stir-bar. In the original set-up a
large circular beaker was used for the outer chamber so that- the
stir-bar was positioned directly below the center of the inner
chamber. It was found that using a set-up such as this would induce a
hydrostatic pressure difference across the membrane due to the vortex?
ing action created by the spinning of the stir-bar. This resulted in
a pressure-driven flow leaving the inner chamber. This could be seen
in that the membrane would bcw outward and the level in the inner
chamber would drop significantly during an experiment. By moving the

stir-bar out from under the inner chamber this problem was avoided.

With both chambers in place distilled water was added to the
outer chamber until the level just reached the level of the membrane.
At this time a step-watch was started and the remainder of the 1600ml
of distilled water and the 400ml of the glycerol solution were quickly
but carefully added to their respective chambers. The two solutions
were added simultaneously at rates such that the fluid levels in the
two chambers remained equal. Stirring in the outer chamber was then
initiated. A Craftsman router with a plastic stir rod powered through
a waerstat variable autctransfcrmer was suspended over the inner
chamber using a second ringstand. The height of the router was set so
that the cross bar at the lower end of the stir rod was about 1/2 inch

above the membrane. Enough power was then applied to the router to

97

provide slow but thorough mixing of the inner solution.

A small sample (about 3ml) of the original glycerol solution was
retained at the beginning of each experiment. Additional samples of
the inner solution were taken at ten minute intervals during the
course of an experiment. A typical experiment would have a total
duration of one hour. Each of the samples was tested for concentra-
tion using an automatic csmcmeter. The csmcmeter is normally switched
on one hour prior to use and calibrated using 100 and $00 milli-csmole

standard solutions as described in the csmcmeter manual.

These samples provided a time history of the concentration of the
inner solution. Ligon. [19]. used the following expression to charac-

terise the concentration of the inner solution as a function of time:

aim - [vi/(viwonoim) +

[Vol (V1+Vo)]c1(0)exp[-tA(D/h) (Hwy/(viva): (4.2 .1)

where 01(t) is the concentration of the inner solution at time t. V
represents volume. A.the surface area of the dialysis membrane. D is
the diffusivity of the solute (glycerol) in the membrane. and h is the
thickness of the membrane. The subscript (1) 'implies the inner

chamber and the subscript (c) the outer chamber.

This expression can be rearranged and the natural leg of each

side taken to yield the expression:

98
1n[[(Vi+Vo)/V°](ci(t)/ci(0)) - (vi/vo)] = st (4.2.2)

where:

s = -A(D/h)(Vi+Vo)/(Vivo) (4.2.3)

For the experimental conditions used in the present work.Vi=4oom1 and

Vo-lGOOml so that:

ln[5ci(t)/4c1(0) - 1/4] = st (4.2.4)

and:
s = 3.125s10” A(D/h) (4.2.5)

Equation (4.2.4) expresses a linear relationship between the natural
log of a simple expression involving the ratio of the concentration at
a given time. t. to the initial concentration and time. Figure 4.2.2
shows the results obtained in a typical experiment conducted as a part
of the present work. These data clearly show that this type of linear
relationship does indeed exist. A linear least squares minimization
was utilized to find the best fit line through this data set. The
resulting slepe is also shown on Figure 4.2.2. The line was not
forced to pass through the origin as would be the case using equation
(4.2.4) exactly. since the start time of each experiment is somewhat
hard to define exactly. This is due to the fact that start-up of an

experiment is rather clumsy and errors of up to a minute in the start

 

3 8
O .0
(I) <1.

0

- ln(

 

 

 

Time (min.)

Figure 4.2.2- Typical results of dialysis membrane
permeability experiment for Cuprophan 100pm
membrane and solute glycerol.

99

100

time are possible.

The quantity (D/h) is taken to be the permeability of the mem-

brane. Solving (4.2.5) for (D/h) yields:

ad = D/h = s/(3.125:10"A) (4.2.6)

where “d will be used to represent the solute permeability of the
dialysis membrane. Using the inside diameter of the inner chamber to

calculate A and inserting this value into (4.2.6) gives:

where if s is in units of (1/min) ”d will have units of (cm/min).
Using the value of s from the experiment shown in Figure 4.2.2 gives

the result:

md = 2.835 x 10" cm/min (4.2.8)

This is compared to the manufacturer's suggested value (see Figure
4.1.3) of (mas-5.2x10-a cm/min. This relatively large difference may be
accounted for by the difference in the temperature at which the two
values were obtained. The manufacturer's value was reported for 37°C
while that generated in the present work was for room temperature (ap-
proximately 23°C). lbdifications to the system which would allow one
to run experiments at various temperatures (both above and below roan

temperature) are being considered for future investigations. These

.‘J

101
modifications will make it possible to more closely evaluate the

correlation between the values obtained using this method and those

reported by the manufacturer.

102

4.3 22251522122112.1322;

This section describes the techniques used while working in the
Bio-engineering Transport Processes (BTP) Laboratory of Michigan State
University. The tasks described are the formation of unilamellar
liposomes ("artificial" cells). and the generation of a series of pho-

tomicrographs using the diffusion chamber.

The liposomes used in the present study were formed from
Leo-lecithin produced by Leon Laboratories of St. Louis. lb. (lot
nnmber 102112). lhen not in use the lecithin was kept frozen in a
plastic jar. About 10 minutes prior to its use the lecithin was
removed from the freezer and placed in a vacuum chamber with dessicant
beads in the bottom and allowed to come to room temperature in this
dry environ-put. Nb v;cuum was drawn on the chamber at this time. A
small sample of the lecithin (about a 1/8 inch diameter ball) was
removed from the bottle and transferred to a glass coverslide. The
ccverslide and lecithin were placed in a 150ml beaker and about 10ml
of a 1:2 (V:V) chloroform-to-methanol solution was added. This mix-
ture was then agitated until all of the lecithin dissolved. The glass
cover slide was then removed from the beaker with tweezers and the
beaker was placed in the vacuum chamber. A.vacuum pump was connected
to the chamber and activated. The pump was left on for 30-45 minutes
causing the chloroform and methanol to evaporate. This process leaves

a thin coating of lecithin on the bottom of the beaker.

The beaker was then removed from the chamber and about 20ml of

103
the solution in which the liposomes are to be formed was carefully
added. Most of the liposomes used in the present study were formed in
a 0.2 mole/liter solution of sucrose. Almost immediately after intro-
duction of the sucrose solution a cloud will begin to form in the
solution. Care was taken not to disturb the solution once this cloud
begins to form. The beaker was covered with Parafilm and placed in a
constant temperature bath preheated to 60°C. The bath was turned off
just prior to placing the beaker into the bath in order to minimize
the disturbance to the solution in the beaker which could be caused by
the turbulent mixing of the bath. The cover was placed on the bath

and the solution allowed to sit for 24 hours.

After the waiting period the beaker was removed from the bath.
At this time one will typically observe that a cloudy ring has formed
in the solution. Experience has shown that the best. most useable
liposomes will be found in and around this cloudy ring. A pipet was
used to draw small samples of the solution (about 5-7ml each) from the
beaker. These samples were placed in centrifuge tubes and centrifuged
at 15.600 G. 15.000 RPM for 15 minutes. After centrifuging the cloudy
material will tend to collect near the top of the tube. Experimental
samples were generally taken from the edges of this cloudy material.
One faces a trade-off as the best liposomes are generally found within
the cloud and yet when sampling from the cloud a significant amount of
"junk" is also obtained. No technique has yet been developed for

separating the good liposanes fran the "junk".

Once the liposcnes (or cells) have been prepared one is ready to

104

generate a series of photomicrographs using the diffusion chamber.
The first step in this process is to prepare the chamber itself. One
must begin by soaking a small (about 2 inches square) piece of the
dialysis membrane in a solution the same as that in which the lipo-
somes (or cells) are in at the onset of an experiment. The dialysis
membrane is soaked for 20-30 minutes prior to use in order to insure
full saturation. A small sample of the liposome suspension is removed
from one of the centrifuge tubes with a micro-pipet and placed in the
sample region of the diffusion chamber with the tap fitting in an
inverted position (see Figures 4.1.1 and 4.1.2). A piece .of .the
dialysis membrane is removed from the solution in which it has been
soaking and stretched firmly by holding at each corner. The membrane
was carefully stretched across the membrane retaining ring such that
no wrinkles are left in the center region. The membrane was released
from the users grasp and instead the corners of the retaining ring are
used to hold both the membrane and the retatining ring itself. The
retaining ring with the membrane clinging to it was then pressed care-
fully over the inverted top fitting such that the sample is disturbed
as little as possible. By holding the retaining ring only. and not
the membrane. the membrane is able to conform to the shape of the top
fitting by slipping between the retaining ring and the sides of the
top fitting. One must be careful to insure that no air bubbles remain
between the membrane and the top fitting. The bulk flow channel in
the diffusion chamber is filled with the initial solution (taken. fees
the solution in which the membrane was soaking) and the tap fitting

pressed into place.

105

The entire chamber was then transferred to the microscope (a
Zeiss Universal Research licroscope) and attached to an X?! traversing
mechanism. One 1/8 inch inside diameter Tygon tube is connected to
the bulk flow outlet port at one end and the other end of the tube was
inserted into a discharge vessel. Another tube is connected to a 35cc
syringe and the syringe and tube are loaded with the new solution to
be introduced into the chamber. All air bubbles are removed frua both
the syringe body and the tube. The free end of the tube is then con-

nected to the bulk flow inlet port.

Provided that there are suitable liposomes in the sample used.
one is now ready to generate a series of photomicrographs documenting
the response of an individual liposome to an induced osmotic imbal-
ance. A Chinon LED Prcmaster 35mm camera with an automatic winder was
used to generate the photos used in the present study. ASA 125 black
and white film was used with an exposure time of 1/30 of a second.
The microscope illuminator was a 60' tungsten bulb with a voltage

input of approximately 12V.

When the cperator is ready to begin the stop-watch is started and
the syringe plunger slowly depressed. Typically the syringe would be
loaded with about 20cc of the new solution and the entire 20cc intro-
duced in about 20-30 seconds. Once the new solution has been
introduced into the bulk flow channel a diffusion process begins
between the solution in the bulk flow region and that on the other
side of the dialysis membrane in the sample region. This diffusion

process creates an osmotic imbalance for the cells in the sample

106

region to which they will respond. Photos are then taken at predeter-
mined intervals. Any special circumstances or observations are noted

in a lab log book for future reference.

It was found to be easiest to use two people to generate a series
of photographs. One person would monitor the chamber through the
microscope keeping the liposome (or cell) of interest in view and in
focus while the other would take the photos at predetermined inter-

vals.

An effective technique was also developed by which a series of
photos could be generated by a single individual. A computer program
written by Tom Gielda. a member of the BTP lab group. when run on the
DEC PUP 11/03 mini-computer in the ET? lab would cause the terminal
bell to ring once every second. Using this program one could count
time by the bells and at the same time keep the liposome of interest
in view and in focus. A stop-watch was also used and checked periodi-
cally to insure that one did not lose count. Typically for the first
1.5-2 minutes a photo would be taken every 10 seconds so that the
operator merely counts to ten and shoots. The time between shots was
then increased to 20 seconds and after 3-4 minutes into the experiment
extended to 40 seconds. A full series of photos would typically have

15-20 individual images taken over a total period of 6-8 minutes.

The film was then processed using standard procedures outlined in
the documentation which comes with the film. Good results were

obtained by placing 8 images on a single 8x10 print using a masking

107
kit. These photos were processed to yield volume as a function of
time using the techniques described in Chapter 4 Section 4 of the
present work. Figure 4.3.1 shows some typical photos of embryos
undergoing a shrink-swell binary flow process with NaCl as the
impermeable solute and glycerol as the permeable solute. Figure 4.3.2
shows some typical photos of liposomes generated using the diffusion

chamber.

The facilities are also available in the lab to record the entire
process on video-tape including a character generator which displays
elapsed time directly onto the tape. The same procedure would be fol-
lowed except that one need not worry about taking photos at specific
times. The character generator would be initialized and set running
.at the beginning of the introduction of the new solution so that no

stop-watch was required. This was the method used by Helkerson in his

investigation of unilamellar liposome permeability. [26].

 

Figure 4.3.1 - Typical ovun photos from diffusion chamber

108

 

Figure 4.3.2 - Typical liposane photos fran diffusion

109

4.4 Imagg Egocgsging Tgchgigues

Cells or liposome ("artificial" cells) were subjected to a
pseudo-step change in extracellular concentration using the diffusion
chamber described in the previous section. The response of the cell
of interest to the non-equilibrium conditions created was documented
through a series of photomicrographs taken at known intervals. These
photo images were processed to yield the volume of the cell as a func-
tion of time. This volume information was then matched to the model
of interest in order to estimate the values of the transport preper-
ties. There are many possible ways to perform the image processing

task.

Shabana. [20]. used slides. rather than prints. and projected the
image of each cell onto a piece of thin tracing paper. The outline of
each cell image was traced onto a separate sheet of paper. These out-
lines were then cut out of the paper and weighed. The ratio of the
weights of each image to that of the first image (at time-0) was taken
as the ratio of the area of the projected images. It was then assumed
that the cell remained relatively spherical so that the radius could
be calculated from the area and from the radius a volume was calculat-
ed. This method assumes that the density of the paper is constant.
It is also very time-consuming as each image must be processed indivi-

dually and by hand.

Another method developed as a part of the present work was to use

a photo-enlarger to project the image from a black and white film

110

111

negative onto a large grid. The diameter of the cell image was then
measured by hand in grid units several times for each image. The
scale of the image could be calculated by measuring the distance
between calibrated scale marks imprinted on each image at the same
time the photo was taken. This method worked fairly well when the
image on the negative was very distinct and the cell remained circu-
lar. Many cases were found in which the image produced by projection
of the negative was not distinct enough to clearly define the boundary
of the cell. This was particularly true for liposome (artificial
cell) images very early in the sequence and for those very late in the
sequence. During the majority of the experiment the difference in
concentration inside and outside the liposome was large enough to pro-
duce a phase-contrast halo around the cell as viewed in the
phase-contrast microscOpe which made the outline'of the liposome eafi-
ly identifiable. However. when. the concentration was nearly in
balance. as in the very early and very late times. this contrast was
not present and the boundary of the cell was not easy to identify in

the negative image.

This problem could be avoided by utilizing positive prints rather
than the negatives. By manipulating the exposure time of the print
the boundary of the cell or liposome could be made distinguishable in
most cases. Unfortunately the images could not be made large enough
to measure accurately without losing resolution in the photo. Instead
the tools of the Computer Image Analysis Laboratory of Michigan State

University administered by Professor Richard Dubes were utilized.

112

The photomicrOgraphic images (typically 8 images to a single 8x10
glossy) were projected onto a video terminal display screen. The pro-
jected image could be made sufficiently large to fill the screen
without losing significant resolution. A. movable cursor was then
positioned at several (typically 20) positions around the boundary of
the cell and the X-Y coordinates of the cursor at each point were
determined by the computer. All of the points entered will have first

quadrant coordinates (positive X and positive Y).

An algorithm called CIRCLE (see Appendix E) was then used to per-
form an integration in radial coordinates to determine the projected
area of the image as defined by the 20 input points. The equivalent
radius of the image was then calculated as the radius of a perfect

circle with the same area as that determined for the image.

The first step in this procedure was to change to an XrY coordi-
nate system whose origin is within the confines of the point set. The
new origin is arbitrarily placed at the location (in terms of the ori-

ginal coordinate system):

X. = (Xmax + xminHZ'o
(4.4.1)

Y. = (1..., + Imp/2.0

where Xm‘x and 1": are the largest x and T values respectively
contained in the data set. and Skin and 1min are the smallest such

values. Thus if the new coordinates are designated 1' and Y' then:

113
x' - x - x,
(4.4.2)

Y'-Y-Y.

These coordinates are then transformed into radial coordinates such
that:
r3 . 1" + Y" (4.4.3)
0 - arctan( X'IY' ) (4.4.4)

For a planar are the area bounded between the origin and the arc can

be expressed as:

I 8 ' ' e
A O‘Io .Ir r dr d9 (4.4.5)

A - 0 Je’ (1/2) :‘49 (4.4.6)

1

In this algorithm r is assumed to be a quadratic function of 9 such

that:
r- e+b0+cO’ (4.4.7)

Substituting (4.4.7) into (4.4.6) and performing the integration on

theta yields the expression:

A - (1/2)[.'o + 289‘ + (2.o+b‘)e’/3 + ch‘l4 + c'e‘ls 9 ]°’ (4.4.8)

1

114
Successive groups of three consecutive points each are used to evalu-
ate the constants a. b. and c and the value of the expression (4.4.8)
is evaluated for each are generated. In practice the (radial) coordi-
nate axes are rotated prior to the calculation of a. b. and e such
that O'-0 for the first point in each group. Thus the value of the
lower limit in the integration on theta is always zero and equation
(4.4.8) need only be evaluated at one value of theta for each point
group (that of the new angle after rotation of the third point in the
current group). The contributions to the area calculated for each are
generated from a group of three points are summed to yield the total
area. In practice each successive pair of points is used in two arcs
generated from the groups of three points (see Figure 4.4.1). Because
of this overlapping of arcs the contribution for each arc is halved
effectively averaging the contributions. Thus if 20 data points are
entered then 20 arcs generated from.three successive points each are

used to calculate the total area of the image.

Thus the user must enter the data points at relatively equally
spaced intervals around the circle. in a clockwise rotation. and in

sequential order around the border of the image.

The scale of the image was determined by positioning the cursor
at two points on the calibrated scale appearing in each image. This
is done prior to entry of the points around the boundary of the image
when using the routine TAIEPT (see Appendix E). The distance between
these two points was then calculated and output with the data points.

By developing all of the images at the same scale with the

 

JZJpohus.mxzenunedammlan
amidtmuy'mmnmer"cmo&mn

 

Saxmd:nxfion:hfingnflxfl. Emma nxfion1hu£gnflxflt

Figure 4.4.1- Point selection sequence and integration
region selection sequence.

115

116

photo-enlarger and not changing the focal distance of the projected
immges the scales determined for each image should ideally be the same
so that each of the individual scales can be averaged to produce the

final scale for all images.

This process was acomplished through two sub-programs. The first
is called TAKEPT and was used to take the data points from the image.
The second is called CIRCLE and does the actual fitting of the points
to the circles. These routines are presented in Appendix E of the
present work. The output from TAIEPT consists of the two scale end-
points. the number of data points entered around the edge of the image
(currently 20 by default). and the actual 1+! data points. This out-
put is written into a file called POINTS.DAT. The routine CIRCLE will
read all the values from this file. perform the integration to deter-

mine the area. and then output the area calculated. the equivalent

radius. and the scale length to a file named POINTS.OUT.

This procedure seemed to work well for most cases. The results
for a given image were reproducible to within 1 percent for a clear
image and if the points were entered carefully. The times required
for image processing using this method were much shorter than those
required for hand processing. By processing each image 2-3 times an
average value could be determined for the radius of each image and
individual errors minimized. The fitting routine. CIRCLE. was very
fast requiring about as much time to execute as it takes to enter the
next command to the computer. Thus a series of 20 images would typi-

cally take about one hour to process using this procedure.

5.1 Ogotic Shrink!“ g; M Embgog

The results reported here were generated using the program MARBOX
to process data reported by Shabana. [20]. on the osmotic shrinkage of
hamster embryos. These results serve to demonstrate the workability
of IlARBOX in determining the permeability of a membrane to water using
the simplified Kedem and Katchal sky model. The raw data reported by
Shabana are tabulated in Appendix F of the present work. The data are
reported in this form in order to provide a set of reference values
for use by other investigators in testing of other computer routines.
In the course of this work one of the problems encountered was the
lack of tabulated data' available with which to test the programs

deve10ped.

Shabana reported numerical values for the nondimensional volume.
V.=V/V.. where V is the volume of the cell as a function of time and

V. is the initial cell volume. of an individual cell as a function of

time. Four mabryos where tested in separate experiments. These
embryos were made to undergo an osmotic shrinkage using the diffusion
chamber described in Chapter 4. Section 1 of the present work. A step
change in the concentration of the solute sodium chloride. to which
the cells‘ are impermeable. was introduced in the bulk flow region of
the chamber and diffused into the sample region through an Enka
dialysis membrane. type 80pm.' All of the experiments were conducted

a t room temperature.

Shabana processed the data by hand using a closed-form solution

117

118

to the simplified Iedem and Katchalsky equation for this case
developed by Terwilliger and Solomon. [21]. This solution assumes
that the cells experience a step change in extracellular solute con-
centration at time equal to zero. Iith the diffusion chamber this
assumption does not reflect the true conditions as the dialysis mem-
brane. which separates the bulk flow region from. the sample region.
will cause the step-change induced in the bulk flow region to be gra-
dually introduced intc the sample region. Tbrwilliger and Solomon's
solution also assumes that the surface area of the membrane remains

constant throughout the experiment and is set to the initial value.

The effects of these two assumptions on the permeability values
calculated were explored using MARBOK and the data for cell 1. The
program IMRBOX has the ability to simulate the conditions assumed in
Terwilliger and Soloaon's solution. It can also take into considera-
tion the lag introduced by the presence of the dialysis membrane in
the diffusion chamber. The concentration of the sample region as a
function of time is calculated using the permeability of the dialysis
membrane. “d4 to the solute present and the algorithm derived in
Chapter 4. Section 1 of the present work. The routine IAIBOX can also
be made to calculate the surface area of the membrane as a function of
the volume of the cell assuming a spherical shape for the cell. These

results are summarised in Table 5.1.1.

Note that the values for the permeability of the cell membrane to

water. LP’ are reported in units of (place). The conversion to these

units from.uuits typical of the ledem and Iatchalsky definition of

*
Table 5.1.1 - Results for Shabana's cell #1.

*

(I)

d

 

Final sum of

 

 

 

 

 

 

Run # conditions L
P squares
l Shabana's value 16.17 -
2 Area constant _3
_ 16.69 8.93 x 10
0d - 1000.
3 Area constant . _3
"d = 1.167 x 10-3 25.69 5.63 x 10 .
4 Area = function(V) _3
_ 18.72 10.81 x 10
o — 1000.
d
Area = function(V) _3
5 3 28.51 5448 x 10

 

«d = 1.167 x 10'

 

 

 

 

reported in units of cm/sec, Lp

119

in microns/sec.

120
solute permeability. (cm’ldyne-sec). (see equation 2.3.24) is made by

using the expression:

Lp(u/sec) = LP(cm'/dyne-sec) [ RT/v' ] x 10‘ (5.1.1)

The first value reported in Table 5.1.1 is the value calculated
by Shabana. The second value was calculated using MARBOX and assuming
that the surface area remained constant and by setting the dialysis
membrane permeability. “d' to a value of 1000(cm/s). Using the
present version of the computer routine which calculates the concen-
tration of the sample region if the dialysis membrane permeability is
set to a value higher than 998.0 (cm/min) the sample region concentra-
tion for all times is set to the value of the concentration of the new
bulk-flow solution. This simulates a step-change in the concentration
of the sample region. The actual value of the dialysis membrane per-
meability is typically 10’1 to 10" (cm/min) so that it is safe to
assume that when the user inputs a value as large as 10’ that they
wish to simulate a step-change response. These conditions match those
assumed in the solution used by Shabana. Cauparison of these two
values show only a.3$ difference between. the calculated permeabili-
ties. This slight difference can be attributed to the difference in
the weighting of the data applied in the two methods of minimization
and to numerical evaluation of the Redem and Ratchalsky equations in
IARBOE. In Terwilliger and Solomonfs solution the data is linearized
and a least squares linear regression used to fit a straight line to

the data. In IARBOK.a least squares is performed directly on the data

121

with no linearization. The closeness of these results support the

conclusion the MARBOX is indeed executing properly.

The third value in Table 5.1.1 differs from the second value only
in that the third value was generated using the manufacturer's
recomended permeability value for the dialysis membrane of
wd=l.1667xlO-'(cm/sec). The surface area was again assmned to rmnain
constant. The results show a 545 increase in the permeability value
calculated for the cell membrane. This result is consistent with the
expected results as if one assumes a step-change in extracellular con-
centration then. one will be consistently over-predicting the actual
extracellular concentration due to the lag induced by the dialysis
membrane. This will force the routine to under-predict the cell mem-
brane permeability in order to maintain the same flux rate at a higher
concentration difference. The magnitude of the change was somewhat
surprising. however. as initial investigations on the diffusion
chamber seemed to indicate that the permeability of the dialysis mem-
brane to the solute sodium chloride was high enough so as to have
little effect on the response of the cells. These results show that a
small effect on the cell volume response can induce a large change in
the calculated permeability values. Thus accurate characterization of
the dialysis membrane permeability to various solutes will be vital to

the future successful use of the diffusion chamber system.

The final values in Table 5.1.1 include the calculation of the
surface area of the membrane as a function of cell volume assuming

sperical geometry. The fourth value used the step-change simulation

122

of the sample region concentration. and the fifth value used the
manufacturer's recommended value for the dialysis membrane permeabili-
ty. These results are consistent with the observation that by
assuming constant surface area one will consistently over-predict the
actual area and hence will cause an under-prediction of the cell mem-

brane permeability.

It is interesting to note that the final predicted values of the
nondimensional volume over time for run 3 and for run 5 were virtually
identical. These two runs both used the manufacturer's permeability
value for the dialysis membrane but for run 3 the surface area was
assumed to remain constant and for run 5 the surface area was calcu-
lated as a function of volume. While they permeability values
calculated differed significantly the final fit of the model to the
dita was not significantly different. This would imply that one
should maintain consistency when calculating predicted responses from
permeability values. That is if the permeability values were calcu-
lated assuming constant surface area then the predicted responses
should also be calculated assuming constant surface area. If the per-
meability is calculated assuming a varying surface area then so should
the predicted responses. One should be aware of this effect when
utilizing published permeability values in simulating cell responses.
As long as one remains consistent with the method used to generate the

permeability values a good fit should result.

The sum of the squares values reported in Table 5.1.1 can be used

as a relative measure of the closeness of the fit between the data

123

points and the predicted curve. A lower sum of the squares value
implies a better fit. From these values it is apparent that the per-
meability of the dialysis membrane will have a major effect on the fit
of the model to the data. The best fit for cell 1 was obtained in run
5 which has the cell surface area being calculated as a function of
the volume and uses the manufacturer's value for the dialysis membrane
permeability. The predicted response for runs 4 and 5 are presented
in Figure 5.1.1. This figure illustrates the effect of the dialysis

membrane on the predicted response.

In an attempt to check the value of the dialysis membrane permea-
bility recommended by the manufacturer the modeling subroutine in
lARBOK.was modified to allow the dialysis membrane permeability to
float as a second parameter. Each of the four data sets reported by
Shabana were processed using this modified routine. The results are
presented in Table 5.1.2. Note that for cell 4 the program was unable
to meet the convergence criterion after 60 iterations. This implies
that this data set does not have a distinct minimum for this model.
The parameter search was varying the parameter values in the third
significant digit around the values reported. One can assume that the
minimum exists somewhere in the near neighborhood of these values if

it exists at all.

The final values calculated for ”d: the dialysis membrane permea-
bility. were reasonably close to the manufacturer's value and fell to
either side of that value. This would tend to support the validity of

the manufacturer's values for permeability as well as the sample

 

1.:Hsv mafia

 

 

 

 

ohm _ cum o.H o.o
_ _

— )— CDC

(7
(ea

0
o no 0 no no
C I>.o
O
o o
o

monocolmmum mcflfismmm

mcounEoE mflmwaofip mcwuopflmcoo

 

 

*
V

Figure 5.1.1- Measured and predicted response of cell #1.

124

Table 5.1.2 - Results letting "d float

at
as parameter.

 

 

 

 

 

 

 

“d = 1.167 x 10.3 "d as parameter #2
“d = 1.247 x 10-3
cell #1 Lp = 28.51 Lp = 27.63
(ed = 0.770 x 10-3
cell #3 L = 39.13
P L = 61.52
9
_ -3
cell #4 L = 27.91 “a " 0'97 x 1°
P Lp = 31.81

 

* 0 reported in units

of cm/sec. Lp

in microns/sec.

Table 5.1.3 - Final permeability values as determined

**

by parameter estimation routine.

k

 

 

 

 

 

 

 

Value reported alue calculated
by Shabana using MARBOX
cell #1 16.39 28.51
cell #2 20.00 46.72
cell #3 19.17 39.13
cell #4 16.17 27.91

 

 

 

** A11 L values reported in units of microns/sec.

12S

 

126

region concentration algorithm derived and used as a part of the
present work. The effect on the calculated permeability of the cell
membrane was very dramatic. This demonstrates again the need to care-
fully characterize the dialysis membrane permeability for the solutes

of interest.

The final recommended permeability values for each of the four
cells as determined by the parameter estimation routine are presented
in Table 5.1.3. These values were generated using the manufacturer's
recomended permeability values for the dialysis membrane and by calcu-
lating the surface area of the cell as a function of cell volume
assuming a spherical geometry for the cell. It is interesting to note
the fairly wide variation in the permeability values calculated.
These cells appear to be identical visually and yet their permeability
to water varies greatly. This points out the advantage of the diffu-
sion chamber in the investigation of population distribution
information through the observation of individuals within the popula-
tion. The final simulated response for cells 2. 3. and 4 are shown in

Figures 5 .l .2 through 5 .1 .4 .

l
2.0

Time (min.)

 

 

 

I
0.0

£3
E\

O

0.94—
8

I
o
e
H

l
(x
(6

0.0

v = V/vo

Figure 5.1.2- Measured and predicted response of cell #2.

127

 

O
o
o
o
l I l l I
E 1 1 l l
.4 :5 <5 6 o’
* .—
V — V/Vo

Figure 5.1.3- Measured and predicted response of cell #3.

128

 

Time (min.)

 

 

 

 

n
D
O
c --
c3
«4
E
O
_- a)
H E
«4
a
'O
l l I I r* 1%?"5
O 0‘ m l‘ ‘0 In
.-( o' o' o' o' :5
*
V - V/VO

Figure 5.1.4- Measured and predicted response of cell #4.

129

5.2.§;g;;y 21g; 1; Ugfertilized Eggster‘ggg

Photomicrographic images were taken of two unfertilized hamster
ova undergoing a process of binary flow. Each ovum was taken from
separate sample groups and tested in separate experiments using the
diffusion chamber. Both of the experiments were conducted at room
temperature and within 3 hours of rmnoval of the ova from the host
hamster. The ovum designated 5 was subjected to a change in extra-
cellular concentration of the permeable solute glycerol from 0.0 to
0.20 (osmol/kg). Ovum 6 was subjected to a change in glycerol conr
centration from 0.0 to 0.25 (osmol/kg). Both ova were initially in a

0.3 (camel/kg) saline solution.

The resulting photomicrographic images were processed using the
computer image analysis techniques described in Chapter 4 Section 4 of
the present work to yield the volume of each ovum as a function of
time. These data were then used to test the execution of the parame-
ter estimation routine for the three parameter case. The parameter
routine was able to converge for both sets of data. Figure 5.2.1
shows the measured and "best fit" predicted response for these two

0V8.

Two additional ova were tested by Th. [29]. and recorded on
video-tape. Both ova were subjected to a change in glycerol concen-
tration from 0.0 to 0.25 (camel/kg). Diameters of the ova were
determined by direct measurement of the projected image from a televi-

sion screen. One of these data sets was successfully processed using

130

 

 

 

D . (—
O O
.._s
\D
I -
I,
u?
I :2 “-
o 5
>
D 0 -—°.
m'
0 O
‘3 "'7:
In
at: . U
5 0' m
0 a D )F—o0
N
[3
0 II
D I.
O [3
O [3 o
O [3 L1;
0 [3
(3 (3
v [:I '-
. =
y:
a _‘3
I T I j I 1 I I I f l &‘ o
1.0 0.9 0.8 0.7 0.6 0.7
*
V = V/Vo

Time (min.)

Figure 5.2.1 - Measured and predicted response of ova #5 and #6.

131

132

the parameter estimation routine. Figure 5.2.2 shows the measured and

"best fit” predicted responses for that ovum. designated ovum 8.

In the processing of each data set using the parameter estimation
routine the parameter search would consistently drive the value of the
hydraulic permeability to the upper bound set for that parameter.
Thus the significance of the calculated hydraulic permeability values
must be questioned. For this case it is believed that the dialysis
membrane itself is rate limiting on the system. due to its relatively
low permeability to glycerol (3.8 microns/sec). and that the ovum
remained in a quasi-equilibrium state. with respect to the water con-
tent. during the initial period of rapid shrinkage. It should be
noted that this should not effect the validity of the calculated gly-
cerol permeability values as the solute permeability will‘be primarily
determined by the rate at which the cell recovers during the later
stages of the experiment. Since the ova solute permeability. values
are much lower than the solute permeability of the dialysis membrane
the model should produce an accurate prediction of the solute tran-
sport process for all times. It is only the water transport that is
limited by the slow response of the dialysis membrane and an upper
bound to the rate of water tranport is quickly approach for reasonable

hydraulic permeability values (30-40 microns/sec).

The solute permeability values calculated were 4.2x10”,
3.0x10”. and 2.9xlO" (microns/sec) for ova 5. 6. and 8 respec-
tively. Jackowski. et al. [30]. reported glycerol permeability values

for mouse ova at room temperature of 1.7x10” (microns/sec). This

.a'
4.4

1!)

Time (min.)

 

 

 

1.0 0.9 0.8 0.7 0.6
*
V = V/Vo

Figure 5.2.2 - Measured and predicted response of cell #8.

133

134

would indicate that the hamster ova are much more permeable to gly-
cerol than are the mouse ova. This is somewhat surprising as one
would expect these two systems to be quite similar. Further work

should be conducted in order to substantiate these findings.

For both ova 5 and 6 the calculated value of the reflection coef-
ficient. c. was 1.0. For ovum 8 the caculated value of c was 0.76 .
The validity of these values is questionable as the reflection coeffi-
cient characterizes solvent-solute interactions and the hydraulic
permeability. and house true rate of solvent transport. is unknown due

to the rate limiting problem.

One significant difference between the data generated using the
techniques of image analysis and that generated by direct measurement
of cell diameter can be seen quite clearly in Figures 5.2.1 and 5.2.2.
That is that the image analysis technique produces "smoother" data.
The fluctuation in the data obtained by direct measurement is due to
the difficulty in determining the diameter of a cell which deviates
even slightly from a spherical shape and to the poor resolution
obtained with use of a ruler. The image analysis technique performs
an integration to determine the projected area of the image and deter-
mines an equivalent radius for a circle with the same area. This
method thereby accounts for deviations from a spherical shape more
accurately than direct measurement. The resolution of the grid on the
computer display terminal used to enter points around the cell boun-

dary is also much greater (spprcxamately 1 part in 400 for a typical

image). This results in more accurate and "smoother" data.

135
It should be noted that the data obtained using image analysis

were checked by using a photo-enlarger to project several of the photo
images onto a grid and measuring the diameter directly. These results
were consistent with those obtained using image analysis and indicated
that no biasing of the data occured through the image analysis pro-

cedure itself.

The tabulated values of the nondimensional volume of ova 5
through 8 are presented in Appendix B of the present work. Also
included in this appendix are the input conditions for the parameter

estimation routine used to generate the results presented above.

The results presented here confirm that the parameter estimation
routine is capable of handling the three parameter problem. and that
'the image analysis methods yielded accurate data. It will be left to
future investigators to increase the data base and to investigate the

implications of these results.

CHAPTER 6

Conclusions

The results of this study demonstrate the workability of the dif-
fusion chamber system in determining the passive transport properties
of an individual cell membrane. This work has also shown. that the
parameter estimation computer algorithm is capable of handling models
of 1. 2. or 3 parameters. The routine has reached convergence on sev-
eral sets of experimental data. The techniques of computer image
analysis have also been shown to work well on well defined photo
images yielding accurate. high resolution data.

Four separate data sets reported by Shabana. [20]. for individual
unfertilized hamster ova subjected to an osmotic shrinkage process
have been successfully processed using the parameter estimation
routine. One conclusion drawn from these results was that the effects
of the dialysis membrane on the concentration history of the sample
region must be taken into consideration when processing data obtained
from the diffusion chamber. Shabana processed the same data sets by
hand using a closed-form solution to the simplified Redem and Ratchal-
sky equations. This solution assumes that the cells experience a step
change in extracellular concentration. When using the diffusion

chamber this assumption does not reflect the true experimental condi-

136

137
tions. The presence of the dialysis membrane between the bulk flow
region and the sample region of the chamber introduces a lag in the
response of the sample region to the step change in concentration
induced in the bulk flow region. As expected it was found that when
processing the same data sets using the relationship derived in
Chapter 4 Section 1 of the present work to calculate the concentration
of the sample region as a function of time the resulting hydraulic
permeability values calculated were consistently higher than those
calculated using the closed-form solution. For the four data sets
processed Shabana reported hydraulic permeability.values in the range
of 16.2 to 20.0 (microns/sec) using the closed-form.solution. Using
the parameter estimation routine with the concentration algorithm
included the same data sets yielded hydraulic permeability values in

the range of 27.9 to 46.7 (microns/sec).

In order to test the parameter estimation routine on a model of
three parameters and to test the image analysis methods as well binary
flow experiments on unfertilized hamster ova were conducted. Two ova
were subjected to changes in extracellular concentration of the perme-
able solute glycerol and their responses were documented via a series
of photomicrographic images. A. third ovum was tested by another
investigator also using the diffusion chamber. The response of this

third ovum was documented on video-tape.

The photo images were processed using the image analysis techni-
ques described in Chapter 4 Section 4 of the present work. The

accuracy of these data was checked by directly measuring the diameter

138

of several of the photo images and comparing the resulting predictions
of volume at various times. The two methods gave similar results
although the image analysis results were much more consistent and
showed the expected trends with much less fluctuation. From this it
was concluded that the image analysis technique yielded more accurate
and more consistent data than did the direct measurement method. The
time involved in using the image analysis techniques was longer
however due to the need to process the film. The actual measurement
process for the image analysis was approximately the same as the

direct measurement process.

Each of these data sets was- successfully processed using the
parameter estimation routine linked to the Kedem and Iatchalsky model
for coupled binary flow in a membrane. In each case it is believed
that the permeability of the ovum to water was so high that the
dialysis membrane permeability to the solute glycerol became the rate
limiting factor in the diffusion chamber system. This could be
observed in that as the hydraulic permeability value increased above a
value of approximately 30.0 (microns/sec) there was no significant
effect on the predicted curve. Thus the values of the hydraulic per-

meability calculated for these data sets have no significance.

This rate limiting of the solvent transport should not have a
significant effect on the calculated solute permeability values as the
dialysis membrane is significantly more permeable to the solute than
are ova. The solute permeability values calculated will be primarily

determined by the rate of volume increase during the later stages of

139
the experiment. The glycerol permeability values calculated ranged

from 2.9x10"3 to 4.2x10"a (microns/sec) for the three cells processed.

CHAPTER 7

Suggestions for Future Work

The results of the present work show great promise for utiliza-
tion of the diffusion chamber in conjunction with various computer
routines in the study of the permeability characteristics of cell mem-
branes. During the course of this work several points which should be

investigated in the future were revealed.

The first of these is the characterization of the dialysis mem-
brane permeability. In the processing of data obtained using the
diffusion chamber the permeability value of the dialysis membrane to
the solutes used was found to play a key role in the resulting cell
membrane permeabilities calculated as describe in Chapter 5 of the
present work. Thus it will be vital to know accurately the permeabil-
ity of the dialysis membrane to each of the solutes of interest. As
future research plans include modification of the diffusion chamber to
accommodate a temperature regulating system one will need to know the
pemeability characteristics of the dialysis membrane as a function of

temperature as well.

The modification of the diffusion chamber to accommodate a tem-
perature regulating system is. in itself. another project which should

140

141
be undertaken in the future. As the overall research effort in the
BTP lab centers on the investigation of the effects of cryopreserva-
tion procedures on various cell types,investigation of cell transport
properties as a function of temperature becomes a key concern. In
order to accurately predict the response of a particular type of cell
to a particular freezing protocol one will need to know the cell's
permeability characteristics over the entire range involved in the

freezing protocol.

Another project which should be undertaken is the modification of
the bulk flow system. The present system of pressurized bottles has
several drawbacks. One is that as flow is introduced into the diffu-
sion chambers bulk flow channel the initial surge can cause severe
disturbance in the sample region. This can make it difficult to keep
the cell of interest in view and in focus. This problem might be
solved by lowering the pressure in the bottles and regulating it more
closely. This would require that a low pressure regulator be
installed in the air supply line. Another problem with the present
bottle system is that due to the way the bottles are interconnected a
relatively long section of tubing exists between the Y-junction which
connects the two bottles to a single pump and the bulk flow inlet port
of the diffusion chamber. This creates a lag of unknown duration
between the time the pump is activated and the time the new solution
actually enters the bulk flow channel. This makes it rather difficult
to clearly define time equal to zero in an experiment. A Y—junction
closer to the inlet port. and use of a separate pump for each bottle.

would probably resolve this problem. While the technique utilizing a

142

manually depressed syringe does not encounter either of these problems
it is more cumbersome to work with and was found by lelkerson. [26].
to be less reproducible. Thus a modified bottle system is the prefer-

able option to pursue.

During several runs using the diffusion chamber it was found that
a sudden rippling of the dialysis membrane would occasionally cause
the entire contents of the sample region to be flushed out of the sam-
ple region. It was also found that when working with the relatively
small liposomes that they had a tendency. on occasion. to drift about
and sometimes squeeze between the dialysis membrane and the top fit-
ting of the diffusion chamber. A. modification of the membrane
retaining ring would probably solve both of these problems. It might
be useful to replace the retaining ring with a cup-like fitting which
would then sandwich the dialysis membrane firmly between two rigid
supports. A small hole in the center of the cup corresponding to the
hole in the spacer ring (which makes up the sample region) would allow
the bulk flow to come into direct contact with the dialysis membrane.
The base of the cup should be made as thin as possible in order to
prevent the region within the hole from becoming an isolated pocket
bypassed by the bulk flow. As the hole in the center will allow view-
ing of the sample region this cup-retainer could be made from most any
material. opaque or transparent. A retainer such as this would pre-
vent rippling of the dialysis membrane and press the membrane finely

against the top fitting preventing leakage from the sample region.

Iany points regarding the statistical nature of the parameter

143
search exist which should be explored in the future. Tb the best
knowlege of the author very little investigation into these aspects of
the parameter search applied to the problem of passive cell membrane
transport has been performed. Enough work relating to this area
potentially exists to justify devotion of an entire thesis project to
just this problem. For instance it would be interesting to apply a
sequential method of minimization to this problem. That is. a method
which adds data points one at a time and adjusts the parameter values
with each additional point. This type of routine can often provide
interesting insights into both the data being processed and the model
being used to simulate the data. It would also be interesting to map
out the sensitivity coefficients as a function of time for the final
parameter values as this can often provide insights into the model
which are not obvious from inspection of the modeling equations. For
instance Papanek's observation (quoted in Chapter 3 Section 7 of the
present work) that in a binary flow shrink-swell situation the initial
rate of shrinkage will be almost entirely dependent on the value of
solvent permeability. and that the shape of the minimum will be pri-
marily dependent on the value of the interaction coefficient. o. and
that the rate of swelling during the final stages will be dependent on
the solute permeability value. should be reflected in the sensitivity
coefficients. Thus one would expect the sensitivity coefficient
relating to solvent permeability to have relatively high values during
the period of initial shrikage and to drop off thereafter. Similarly

for the other sensitivity coefficients.

The investigation of models other than the Kedem and Katchalsky

144

permeability model used in the present work is an obvious area of
research for the future. This is particularly true with regards to
using the Johnson and Wilson power series solution. outlined in
Chapter 3 Section 6 of the present work. as a generator of starting
estimates to the parameter values. Use of this model should improve
the initial estimates and hence reduce the total computational effort
required. The Papanek model. described in Chapter 3 Section 7 of the
present work. also needs to be investigated. This could be done in
conjunction with the investigation of the statistical nature of the
parameter search applied to the passive transport problem. Various
methods exist by which one can statistically compare two or more
models. It would be interesting and very useful to perform such a
comparison between the Papanek model and the more commonly used K-K
model. It would also be interesting to apply the Kedem and Iatchalsky
resistance model. particularly to the dialysis membrane where one
could assume a capillary structure and thereby simplify the model con-

siderably.

A final area which needs further work is the deve10pment of image
analysis methods. Reasonable results were obtained in the present
work using very simple methods of analysis. lore sephisticated
methods with better reproducibility should be developed as the com-
puter image analysis techniques have the potential for becoming an
elegant and easily utilized method of photo data processing. This

project in itself could become a major undertaking.

The present work has demonstrated the workability of the diffu-

145
sion chamber as an experimental tool. and the computer programs
written as a part of that work as analytical tools. It is now up to

future investigators to refine and expand these techniques.

APPENDICES

146

APPENDIX A

Subprogram Unit HARBOR

The ordinary least squares based parameter estimation routines
main driving program is contained in a file name lARBOX.FOR. This
file also contains three support routines called only from .IAIN.
which will also be described here. These support routines are BOUNDS.

IIINVER. and DETERI.

The input to the routine is through a data file. Input is read
from logical unit 2 so that under the RT11 operating system the input
file has the name FTN2.DAT (see R111 FORTRAN IV Users Guide page 3-5).
All input is in the free format mode so that no special formating is
necessary. Input values on a single line should be seperated by com-
mas. The order of input and the definition of each of the input

variables is included in the program listing.

The main routine embodies a variety of parameter estimation
routines all based on the ordinary least squares (OLS) method. The
form of the routine to be executed is set by the user through) the
input variables BOEFLG. and IFLAG. By setting one of these variables
to 1 the user introduces modifications to the basic OLS routine. If

both of the flags are set to 0 then the routine is the 0L8 method with

4‘

147

a modification which checks for violation of the parameter bounds set
by the user during input. The OLS routine will have the fastest exe-
cution time per iteration on the parameter values but will probably

require more iterations to reach convergence than the other forms.

The first modification allows the user to set upper and lower
limits on the values of each of the parameters. In the current ver-
sion of the program this option is invoked by default. Within the
input section of the program there is a set of statements which can be
easily modified to allow the user to make this option selectable via
the input variables. This is not recomended however as this option
adds very little to the execution time of the program and will prevent
the search from diverting to unrealistic values of the parameters
(such as o<0 or o)1.0). For instance in the Ir! equations negative
values for any of the parameters will result in meaningless solutions.
Thus the user should set a lower limit on each of the parameters of no
less than 0.0 . After each iteration the newly calculated parameters
are compared to the bounds set by the user. If the new value violates
either of the bounds then the value of the parameter is set to the
value of the violated bound and execution continues. The checking is

done in subroutine BOUNDS.

The second option is invoked by setting lFLAG-l. This causes the
routine to execute as Narquardt's method (see Chapter 4 Section 5 of
the present work). larquardt's method modifies both the size and

direction of the step in parameter values taken in each iteration.

148

The final option is invoked by setting BORFLG to one. This
causes the routine to execute under the Box-Kanemasu minimization

method.

Note that the Box-lanemasu method is not compatible with lar-
quardt's method so that one should not set both BOXFLG and IFLAB equal

to one.

The convergence of the routine is based on the change in the
parameter values through the input variable TOLER. Execution is ter-
minated and final values of the cell volume with time calculated when

the condition:

ABi’Bi ( TOLER

is met for all of the parameters (i-l.2......P).

The support routine am is a subroutine which will invert an
(nxn) matrix where (n53). Subroutine DETERN calculates the determi-

nant of an (nxn) matrix where (n53).

In utilizing this routine the user should have patience. The
results in a parameter estimation routine can be unexpected and the
corrective actions required to enable convergence are often learned
only through experience. By utilizing the various options available
the user can run a particular data set under several different pro-
cedures. For instance the user can "trick" the routine into

manipulating only one of the parameters in the three parameter model

149
by setting the upper and lower bounds of the other two parameters to

the value of the initial estimate and invoking the bounding option.
Thus one can "zero in" on the final parameter values one at a time.
This is often helpful when the initial estimates are not very good and
the routine has trouble converging. By allowing only one parameter to
vary at a time the user can often get an idea of the direction in
which to change the initial estimates in order to approach the solu-

tion. Even this method will not work in every case. however.

 

 

 

 

OGOOOOOOOGOOOOOGOOOOOOOOOGOOOOOOOOOOOOOGOOO0000606630GOG

150

PRmRAM REFEREVCE: BECK. JABES V.; ARNQD. KENNETH 1.3
PARADETER ESTIMATION IN ENGINEERING
AND SCIENCE. JOHN WILEY _SONS. INC. .
1977. CALL TA340.B39.

THIS PRMRAM IS A PARADET‘ER ESTIMATION RWTINE BASED ON
BAS- ON THE ORDINARY LEAST SQUARES RWTINE. IT ALSO HAS
AN OPTION FOR CONSTRAINING THE PARABT‘ER SEARCH.

IT IS ALSO SETUP TO RUN AS MARQUARUI'S RWTINE OR AS THE
BOX-KANEMASU IETHOD.

INPUT TI) THIS RWTINE IS THRCIIGH A FILE NABED FINZ .DAT
THIS FILE SHOULD CONTAIN IN THE FOLLWING ORDER:

D-UGJIFLAG.BOXFLG
ACONFL.P

II

T(l) .Y(1)

T(2) .Y(2)

fin) .im)
BETA(l) .moau(1).mmm(1)

BETA(P) . IBBETA(P) .IBBETA(P)
TOLER
DT
CIZERO.CISTEP
CSZERO.CSSTEP
DPERMI.DPERDB.DEPTH
RADIUS
VDEAD
( IF MFLAG=1 (MARQUARDI‘S BETHOD) THE‘I: )
LABBDA.FACLAN

WHERE:

onus INTEGER Dill} WTPUT CONTRG. CHARACTER

BOXFLG INTEGER FLAGS FOR BOX-KANEMASU 1811100
1 = YES
0 = N0

MFLAG INTFBER FLAGS FOR MARQUARUI‘S METHOD
1-YES
O-NO

AmNFL 11(me FLAGS FOR ASSUMPTION 0F CONSTANT DEDBRANE
SURFACE AREA
0 = (INSTANT AREA
1 a Nm-CONSTANT AREA

P INTEGER NUDER OF PARAMETERS IN DIEL

II INTEGER NULBER OF DATA POINTS

T REAL(II) TIDE DATA IN VECTOR FOR]! IN (SEC)
Y REAL(II) DEPENDENT VARIALBLE IN VECTOR FORM

BETA REAL(P) PARADETER ESTIMATES

LBBET‘A REAL(P) [HER BWND ON PARAMETERS

UBBETA REAL(P) UPPER BWND ON PARADET‘ERS

DT REAL AFROXIllATE VALUE OF TILE INGREDIENT FOR

 

(5000000000000000000000000OOOGOOOOOOOOOOOGOO000000600000

151

INTERATION STEP SIZE (SINCE INCREMENT
MAY NOT FIT TOTAL TILE INTERVAL EVELY
ACTUAL STEP SIZE IS CALCULATED 1N "RK4"

CIZERO REAL INITIAL (DNCENTRATION OF IMPERMEABLE SG.UTES
CISTEP REAL NEW IMPERMEABLE SCEUTE CONCENTRATION

CSZERO REAL INITIAL mNCENTRATION OF PERMEABLE SEUTE
CSST'EP REAL NEW PERMEABLE SG..UTE CON(ENTRATION

NOTE- CONCENTRATIONS IN

(MOLES/ CM”3)

DPERMI REAL DIALYSIS MEbBRANE PERMEABILITY TO SEUTE I (CHI SEC)
DPERIB REAL DIALYSIS MEMBRANE TERMEABILITY TO‘ SG..UTE 8 (CN/ SEC)
DEPTH REAL DEPTH OF SAMPLE REIGN IN DIFFUSION CHABBE (CM)
DEPTH =- (TRANSFER AREA) / (SAMPLE VEUBE)
RADIUS REAL CELL RADIUS IN cm
VDEAD REAL FRACTION OF INITIAL vmmm ATRIBUTED TO OSBDTIC
DEAD SPACE. 'IHAT IS SPACE NOT INV(LV- IN
ACTIVE TRANSPORT. ( 0.0 (3 VDEAD => 1.0 )
LADBDA REAL SEE MARQUARDT'S METHOD
FACLAN REAL SEE MARQUARDT' S DBTHOD

D-UG IS A ULTY-LEVE IEBUGGING (DNTRG. GARACTER. ITS VALUE WILL
DETERMINE 'THE LEVE OF D-UGGING EINTWTS GEEATED BY mE EmRAM
DURING ICUTION. A BIG. VALUE RESULTS IN TIME EXTENSIVE PRINTOUTS.
'THIS IS VEY HANDY FOR PRINTING WT INTERMEDIATE VALUES 0F VARIABLES
NOT NORMALLY DESIRED IN 'IHE OUTPUT OR FOR PRINTING A MESSAGE PRIOR
TO A SUBRGJTINE CALL IN ORMR TO DETERMDIE WHERE A PRWRAM IS RUNNING
INTO TRWBLE. A STATEIWT OF THE FOLLWING TYPE IS RECOMENIED:

IF (D-UG .GT. 2) PRINT . . . .

THIS PRUIRAN IS SET UP FOR PR(BLENS OF UP TO 3 PARADETERS
ALTHWGH IT MAY VEY EASILY BE DDDIFIED TO WORK WITH ANY
NUIBE OF PARADE'ERS. TO IN) THIS THE USER MUST CHANGE
'THE DITENSION OF ARRAYS BETA. BSTORE. IBBETA. IBBETA. DELB.

H. X. m. AND CHANGE.
ONE HIST ALSO SUPPLY RWTINE TO CALCULATE ms IETERNINANT AND
INVESE OF THE LARGE MATRIX (PIP) .

THE USE MIST SUPPLY A SUBRGITINE NAIED "DEQBDD" IN 'THE
FOLLWING FORMAT:

SUBRGJTINE IDEQBDD (Y1 . Y2 . T. DYlDT. DY2DT. BETA)
REAL DYlDT. DY2DT. T. Y1 . Y2 . BETA( 3)
( OPTIONAL (DEMANDS: )
INTEGE DEBIK}
COMDDN/ BLOCR9/ VDEAD
CONTON/ BLOClO/ AmNFL
COULDN/ Bm/ DEBIK}

WHERE:

Y1 FIRST DEPENIENT VARIABLE (TYPICALLY V‘)

Y2 SECOND DEPENMVT VARIABLE (TYPICALLY NS‘)
T - REPRESEVTS THE INIEPENIENT VARIABLE (TIDE)
DYlDT DERIVITIVE OF Y1 W.R.T. TIDE

DY2DT DERIVITIVE OF Y2 W. R.T. TIDE

BETA ITINATED VALUE OF PARAIETERS

((‘t'l‘A‘

0001

0002

0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013

0014

0015

0016

0017
0018

OOOOOOOOOOOOOOOOO

0000 000 on

000000

152

VDEAD - REAL VALUE INDICATING FRACTION OF CELL DEVOT'ED TO
OSDDTIC DEAD SPACE ( 0 <= VDEAD ( 1.0 )
AmNFL - INTEGE CONTRG. VARIABLE TO SET AREA CONSTANT
ACONFLIO THE A=CONSTANT
ACONFL-l THEN A=F(V(I.UIE)
D-UG - INTEGE DEBUG CONTRG. VARIABLE (SEE ABOVE)

'IHE VALUES OF Y1. Y2. TIDE. AND BETA ARE INPUT TO THE SUBRETINE
AND SHWLD NOT BE ALTERED.
THE VARIABLES DYlDT AND DY2DT ARE 'IHE SUBRETINE OUTPUT.

PRERANED BY: STEVE NELE
3/ 83

DIDENSION BETA(3) .IBBETA(3) .IBBETA(3) .BSTORE(3)
_.EANGE(3) .DELB(3) .E(50) .ETA(50) .H(3)
_.PR(3 .3) .TIIE(50) .X(50.3) .XTX(3 .3) .XTXM(3 .3)
_.Y( 50)
INTEGE I. I. II. P. !OUNT. L. SFLAG. CDIFF. CFLAG
_. BOXFLG. BFLAG. “FLAG. DIUG. DOWN. A(DNFL
REAL PR. LBBETA. IDDEL. MINVE. LAIBDA
COMIDN/ BLOC!1/ BETA. LBBETA. IBBETA. BSTORE
COMWN/ BLOCKZ/ TIDE. Y. ETA. X
COMIDN/ BLOC!3/DELB
COMIDN/ BIJOC!4/DT ‘ .
CONDDN/ BLOC!5/ RADIUS
CONlDN/ BLOCKGI CSZEO. CSSTEP. CIZEO. CIST‘EP. DPERDB. DPERMI. DEPTH
COMIDN/ BLOCKS/P. CDIFF
COMDDN/ BLOC!9/ VDEAD
CONDN/ BLOClO/ ACONFL
CONWN/ BIB/DEBUG

SET ITEATION LOOP COUNTER
!OUN'T'-0

SET CDIFF TO INITIALLY USE A FOEARD DIFFEREE

CDIFFIO
eeeseeeeeeeeeeeeeeeeeeeeeeeeeeeeee

DATA INNT SECTION:
READ (2.‘)D-UG.|IFLAG.BOXFLG

DEFAULT VALUE EEC!S FOR PARAIETER
BENDS VIOLATION
BFLAG INTEGE FLAGS FOR BENDS EEC! 0N PARADETERS
1 8 EC! FOR VIOLATION
O = NO EEC!
BFLAG-l
READ (2.‘)ACONFL.P

 

 

0019
0020
0021
0022
0023
0024

0025

0026

0027
0028
0029
0030
0031

0032

0034
0035
0036
0037
0038
0039
0040
0041
0042
0043

0044 ,

0045
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059

0000

00

00

00000

10

90

90

90

90

90

92

19

92
17

93

92

927

153

READ (2.")II

DO 10 I=1.II

READ (2.‘)TIBE(I).Y(I)

CONTINUE

READ (2.‘)(BETA(I).IBBETA(I).IBBETA(I).I=1.P)
READ (2.‘)TOLER

SET TOLE FOR USE IN CHECKING FOR NEAR CONVENIENCE
AND SUITE TO CENTRAL DIFF DERIV.

TOLE=10 .‘TOLER

READ IN NUBBE OF STEPS FOR R-! RETINE INTEGRATION
READ (2.‘)DT

READ IN EXPERIIENTAL ENDITIONS
READ (2.‘)CIZERO.CISTEP
READ (2.’)CSZERO.CSSTEP
READ (2.’)DPERMI.DPERDB.DEPTH
READ (2.‘)RADIUS
READ (2.‘)VDEAD

IF MARQUARUT DETHOD READ IN LAIBDA AND DL
IF(MFLAG.E.1)READ(2.‘)LAIBDA.FACLAM

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

DATA VEIFICATTON BLOC!:

PRINT 900
FORMAT(1El. ' SPECIFIED RUN ENDITIONS: ')
PRINT 901.11
FORMAT( ' NUDBE OF DATA POINTS ='.I3)
PRINT 902
FORMATU . ' INTEPENDENT DEPENDENT')
PRINT 903.(TIIE(I).Y(I).I=1.II)
FORMAT(2E15.6)
PRINT ‘. ' INITIAL PARATETER ESTIMATES: '
PRINT 904.(I.BETA(I).I-1.P)
FORMAT(' BETA'.Il.' a '.E15.6)
IF(BFLAG.NE.1) GOTO 19
PRINT 921.(I.BBETA(I) .I.IBBETA(I) .I=1.P)
FORMAT(E15.6.' < BETA'.Il.' ( '.E15.6)
GOTO 17
CONTINUE
PRINT 922
FORMAT( ' NO ENS'TRAINT ON PARAIBTERS')
CONTINUE
PRINT 934.TOLER
FORMATUI.’ MEANE FOR CONVENEE = '.E10.4)
PRINT 926
FORMATU/I . ' EXPERIENTAL ENDITIONS: ' .//)
PRINT 927.CIZEO.CSZERO
FORMAT( ' INITIAL (X)N(ENTRATIONS: ' ./ .
_' INITIAL IMPERIEABLE S(LU'TE ENCENTRATION = '

0

1

2

3

4

1

2

4

6

0060
0061

0062
0064

0065
0067

0068
0069

0070
0072
0073
0075
0076
0077

0078
0079
0080

0081
0082
0083
0084
0085

928 FORRAT(' NEW VALUES OF CONCENTRATION IN RULR S(LUTION: '.I.
_' NFII IMPERMEABLE SILUTR CONCENTRATION = '
_.E16.8.' mLES/CC'J.

_' NEW PERMEABLE smUrF. CONCENTRATION = '
_.E16.8.' IDLES/CC'J)
IF(IIOUNUII.NF.3)FRINT 929.DPERMI.DPERMS.DEPTH

929 FORMAT(' CONDITIONS FOR DIALYSIS MEMBRANE IN DIFFUSION (nAInFR:'
_./.' PERMEABILITY To IMPERMEABLE SCLUTF = '
_.E16.8.' CM/SEC'J.

_' PERMEABILTI‘Y To PERMEABLE smUTF = '
_,El6.8.' CM/SEC'./.
_' DEPTH OF SAMPLE REGION (vaIm/ARRA) = '
_.E16.8.' cum)

IF(IIODNUII.BR.3)PRINT 938

938 FORMAT(' DIALYSIS IFIBRANF FFRIUIADILITY To IIIFFRIUIARLF SILU'IE'J.
_' IIILL FLOAT AS SECOND PARADETER'JI.

_' DEPTH OF SAMPLE RFCION (VCLUlB/AREA) - '.El6.8)
PRINT 930.RADIUS.VDEAD

930 FORMAT(' INITIAL CELL RADIUS = '

_.El6.8.' cum/I.
_' FRACTION OF DEAD SpAdI IN INITIAL vaIa = '
_.F6.4.' (mum/(INITIAL VCLUIm'J)
IF(ACONFL.UO.0)PRINT 931
931 FORMATU.’ MEBBRANE SURFACE AREA ASSUIED CONSTANT')
IF(ACONFL.DO.1)PRINT 932
932 FORMATU.’ ImIBRANR SURFACE AREA IIILL VARY IIITU vaIe')
PRINT 933

933 FORMAT(1E1)

C
COOCOOCCOOOOOOOOOOOOOOOOOOOOOOOOOO
c
c RESIN CALCULATION PROCEDURE:
c
c To START PROCEDURE SET DELB(I) =- 101. (BETA(I))
c THIS IS To SET INCREIENT FOR FIRST EVALUATION OF
c PARTIAL DERIVITIVES IN IoDFL IIFIm ARE BASED ON
c VALUE OF DELB.

Do 18 I=1.P

DELB(I)-0.1‘BETA(I)

18 CONTINUE

c

C .0....#.... TOP OF “IN mmm LOOP OOOOOOOOOOOQ

c
20 CONTINUE
KOUN’D-KOUN'HI
DO 23 I=1.P
RS'mRF(I)-DF.TA(I)

25 CONTINUE

c

c CALL TO USER SUPPLIED SURRoITINF ”IDDEL"

C

154

_JE16.8.' MOLES/CC'.I.
‘_' INITIAL PERMEABLE SOLUTE CONCENTRATION = '
{_.E16.8.' MOLES/CC'./)

PRINT'928.CISTEP.CSSTEP

0086

0087
0088
0089
0090
0091
0092
0093
0094

0095
0096
0097
0098
0099
0100
0101
0102

0103
0104
0105
0106

0107

0108
0109
0110
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128

155

CALL IDDEL( II . TIDE. ETA. X. BETA. UBBETA. LBBETA. 1)
C
C CALCULATE AND STORE SUM OF SQUARES
C
SQUAR#0.0
DO 30 Iil.II
E(I)=Y(I)-ETA(I)
SQUAR=SQUAR+E(I)"2
30 CONTINUE
SSTOREPSQUAR
PRINT'905.SQUAR
905 FORNATU.’ SUM OF SQUARES FUNCTION FOR THESE VALUES ='.E15.6)
C
C BECK'S EQUATION:
C
D0 31 I=1,P
DO 31 J=I.P
m(I.J)=0.0
DO 32 L=1.II
XTX(I.J)=XTX(I.J)+X(L.I)"X(L.J')
32 CONTINUE
XTX(J.I)§XTX(I.I)
31 CONTINUE
C
C SET UP XT'X DUMMY FOR NARQUARUI'S mDIFICATIONS
‘ DO 33 I81,P
D0 33 151.P
XTXN(J.I)=XTX(J.I)
33 CONTINUE ‘
C
C SET COUNTER FOR MAEQUARDTS METHOD
IIOUNTEO

Ott¢¢ttttt0 TOP OF MARQUARDT LOOP OCOOOOOOCOCCO

000

34 CONTINUE
NKOUNT=IKOUNT¥1
IF(IFLAG.NE.1) GOTO 36
DO 35 I31,P
XTXNU. I)=XTX(I. I)‘(1 .0+LADBDA)
35 CONTINUE
36 CONTINUE
CALL HINVER(XTXN.PR.P)
D0 37 I=1.P
H(I)=0.0
DO 38 LF1.II
H‘I)=H(I)+X(L.I)‘E(L)
38 CONTINUE
37 CONTINUE
D0 39 I-1,P
DELB(I)=0.0
D0 40 1.2-1 .P
DELB(I)IDELB(I)+PR(I.L) 'H(L)
40 CONTINUE
39 CONTINUE

 

 

 

0129
0130
0131

0132

0134
0136
0137
0138
0139
0140
0141
0142
0143
0145
0147
0148
0149
0150
0151
0152

0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0165
0166
0168
0169
0170
0171
0172

0173

156

D0 41 I=1.P
BETA(I)=BSTORE(I)+OELB(I)
41 CONTINUE

C CHECK FOR OUT OF BOUNDS

IF(BFLAG.EQ.1) CALL BOUNDS(P)
C
c t O C O O O t O O a O c O O 8 t 8 o C O C t O o t t C O t 8 O O
C

C CHECK FOR BOX-KANEMASU DETHOD AND BYPASS THIS SECTION IF NOT
C
IF(BOKFLG.NE.1)GOTO 50
ALPHA!1.0
A91.1
42 CALL IDDEL(II.TIDE.ETA.X.BETA. UBBETA.IBBETA.0)
SALPHAF0.0
D0 43 I=1.II
SALPHASSALPHA*(Y(I)-ETA(I))“2.0
43 CONTINUE
IF(SALPHAALT.SSTORE)GOTO 46
IF(ALPHA.LE.0.01)GOIO 45

ALPHAFALPHAI2.0
DO 44 I=1,P
BETA(I)iBSTORE(I)+ALPHAfiDELB(I)
44 CONTINUE
GOTO 42
45 CONTINUE
C ERROR CONDITION

PRINT '. ' UNABLE TO REDUCE SUM OF SQUARES IN BOX METHOD'
FRINT".' EHECUTION CONTINUING AT’LAST VALUES OF ALEHA' '
46 CONTINUE
6-0.0
DO 47 I=1,P
DO 47 I81.P
G=G+DELB(J)’XTX(I.J)*DELB(J)
47 CONTINUE
TEST%SSTOREFALPHAFG‘(2-(1/A))
HBOX=AFALPHA
IF(SALPHA. LT. TEST)GOTO 48
HBOX2=G‘(ALPHA“2 .0)/(SSTORE+SALPHA+(2. 0*ALPHAFG))
IF(HBOXZ. LT.BOX)EOX=BOX2

48 CONTINUE
m 49 I81,P

DELB(I)=DOX’DELB(I)

BETA(I)=BSTORE(I)+DELB(I)
49 CONTINUE
C
C END OF BOXPKANENASU METHOD.BLOCK
C
c o t C O o o o O t c t O o t o t o t o t t t O C O o t c t c o t C
C

50 CONTINUE

C
C IARQURDT'S METHOD BLOCK:
C
C AGAIN CHECK FOR MARQUARDT'S METHOD FLAG:

0174
0176
0177
0178
0179
0180
0181
0183
0185
0187
0188
0189

0190
0191
0192

0193
0194
0196
0197
0198

0199

0201
0202
0203

0204
0205
0206
0208

0209
0210
0211

00000

0000000

000

157

IF(MFLAG.NE.1)GOTO 54
CALL M)DEL( II. TIDE. ETA. X. BETA. UBBETA. LBBETA. 0)
SQUAR=0.0
DO 51 I=1.II
SQUAR=SQUAR+(Y(I )-ETA(I) )“2
CONTINUE
IF(SQUAR.LT.SSTORE)GOTO 53
IF(MKOUNT.GT.10)GOTO 52
IF(DEBUG.GE.1)PRINT 998.FACLAM
FORMATU INCREASING LADBDA BY FACTOR 0F '.E10.4)
LABDA=LAIBDA‘FA(I.AN
GOTO 34

51

998

Ottoman”... 30mg 0F MARQUARDT LOOP OOOOOCCCOOOCO
ERROR CONDITION WTPUT BLOCK:

CONTINUE
PRINT 912
FORMATU.’ UNABLE TO REDUCE SUM OF SQUARES WITH LAIBDA'J.
_ ' AFTER 10 ITERATIONS. CONTINUING EXECUTION')

CONTINUE
IF(D-UG.GE.1)PRINT 999.FACLAM
FORMATU REDUCING LAIBDA BY FACTOR 0F '.E10.4)
LAIBDA'BLAJBDA/FACLAN

mNTINUE

52
912
53
999
54
‘umuomuuuo 30mg. 01: ESTIMATION BLOCK ”tuttouuto
BmIN CHECK FOR CONV‘EWCE:
IF CONSTRAINT PRCBLEN (EEC! PARAIETERS FOR GIT 0F BWNDS:
IF(BFLAG.m.1) CALL BWNDSU’)
ITERATION COMPLETE. CHECK FOR CHANGE IN PARABETER VALUES:
D0 60 I-1.P

CHANGE(I)-ABS(DELB(I) )/(ABS(BSTORE(I) )+1E-16)

60 CONTINUE

C
C BE SURE CDIFF .m. 1 SO THAT AT LEAST ONE IT'ERATION US. CENTRAL

C DIFFEREKE APPROXIMATION FOR DERIVITIVES

C

SFLAGa-l

D0 65 I-1,P

IF(CDIFF.NE.1 .OR. CHANGE“) .GT. TILER)SFLAG-0

65 CONTINUE
C
C (EEC! FOR NEAR CONVERGENB AND SWITCH TO CENTRAL DIFF IF TRUE
C

CFLIO-l

DO 66 I-1.P

IF(CHANGE(I) .GT. MER2)CFLAG-0

0213
0214

0216
0217
0218
0219

0220

0222

0223

0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236

0237
023 8
023 9
0240

66

C

158

CONTINUE
IF(CFLAG.m.1)CDIFF=1

C OUTPUT NEW PARALETER VALUES:

C

906

907
C

PRINT 906.KOUN'T

FORMATU.’ NEW PARAIETER VALUES AFTER '.I3.’ ITI'ERATIONS:')
PRINT 907.(1.BETA(I).CHANGE(I).I=1.P)

FORMAT(' BETA'.Il.' =',m5.6.101. '70 CHANGE ='.E15.6)

C CHECK IF CONVERGBWG BET:

C
80
C

C TOLERANG NOT IET.

C

C
500
C

IF(SFLAG) GOTO 500
RE-ITERATE:
GOTO 20

CONTINUE

C DID PRmRAN SNUENE:

C

908

909

70
935

910

911

CALL ADDEIJ 11. TIDE. ETA. X. BETA. IBBETA. IBBETA. 0)
PRINT 908
FORMAT(1m. ' FINAL ESTIMATED PARABETER VALUES: ' .//)
PRINT 9-093(IDBETA(I)DI.11P)
FORMAT“ ' BETA' . 11. ' =' .m5.6)
SQUAR=O .0
DO 70 181.11
SQUAR'flUAR+(Y(1)'ETA(1))”2 .0
CONTINUE
PRINT 935.&UAR
FORMATKII.’ FINAL SUN OF 'mE SQUARES FUNCTION 3 3315.6)
PRINT 910 '

FORMATU/ . ' INIEPENDENT' .81. 'DEPE‘IIENT VARIABLE' ./ .

_4X. 'VARIABLE' .7X. 'IEASURED' .81. 'DDDEL' .Il)

PRINT 911. (TIDEU) .Y(I) .ETA(I) . I=1 . II)
FORMAT(3E15.6)

STOP

MD

 

 

0001
0002
0003
0004
0005
0006
0007
0008
0010
0011
0012
0013
0014
0016
0017
0018
0019
0020
0021
0022
0023

C

159

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C mIS SUBRCIITINE CHECKS FOR VIOLATION OF THE PARADETER BOUNDARIES:

C

11

10

901
902

SUBROUTINE BOUNDS(P)
DUENSION BETA(3) .IBBETA(3) .UBBETA(3) .BSTORE(3) .DELB(3)
REAL LBBEIA
INTEGER P
COMWN/ BLOCKI/ BETA. LBBETA. IBBETA. BSTORE
CONDON/ BLOCK3/DELB
DO 10 I=1 .P
IF(BETA(I).GT.LBBETA(I)) GOTO 11
PRINT 901 . I
DELB(1)=LBBEIA(I)-BSTORE(I)
BETA(I)=LBBETA(I)
GOTO 10
IF(BETA(I).LT.UBBETA(I)) GOTO 10
PRINT 902.1
DELB(I)=UBBETA(I)-BSTORE(I)
BETA(I)=UBBETA(I)
CONTINUE
RETURN
FORMAT( ' DOVER BGJND VIOLATION BY BETA' .11)
FORMAT( ' UPPER BGJND VIOLATION BY BETA' .11)
END

 

 

 

0001
0002
0003
0005
0007
0008
0009
0010
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0023
0024
0026
0027
0028
0029
0030
0031
0032

160

C
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
C
C THIS SUBWTINE WILL CALCULATE THE INVERSE OF A MATRIX "A"
C WHICH IS NXN AND PLACE THE RESULT IN "B" (N ( 4).
C
SUBRCIITINE N1NVER(A.B. N)
IF(N.NE.1)GOTO 5
IF(A(1.1).EQ.0.0)GOTO 6
B(1.1)=1.0/A(1.1)
RETURN
5 DETIFDETERN(A.N)
IF(DETA.NE.0.0) GOTO 10
6 TYPE ‘.' ERROR - SINGULAR MATRIX'
STOP
10 CONTINUE
DO 20 181.N
DO 20 I=1.N
NN-N-l
DO 30 K81,NN
D0 30 L=1.NN
KK=K
IF(K.GE.J) KK=K+1
LLBL
IF(L.GE.I) LL=L+1
COFA(K.L)=A(KK.LL)
30 CONTINUE
OETCAPDETERNICOFA.NN)
B(I.J)=DETCA‘((~1)“(I+J))IDETA
20 CONTINUE
RETURN
END

 

 

—_.___—__—#-

 

161

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
0001 FUNCTION DETERN(A.N)
0002 DIENSION A(3.3)
0003 IF(N.NE.1) GOTO 10
0005 DETERM-A(1.1)
0006 RE'IURN
0007 10 IF(N.NE.2) GOTO 20
0009 DETERDhA(1.1)‘A(2.2)-A(2.1)‘A(1.2)
0010 RE'IURN
0011 20 IF(N.NE.3) GOTO 30
0013 DETERM- A(1.1)‘(A(2.2)‘A(3.3)-A(2.3)"A(3.2))
_ -A(1.2)‘(A(2.1)‘A(3.3)-A(2.3)‘A(3.1))
__ +A(1.3)"(A(2.1)*A(3.2)-A(2.2)‘A(3.1))
0014 RE'IURN -

0015 30 TYPE 9. ' ERROR - NATRIX IN IBTERM CALL LARGER THAN 3X3'
0016 STOP
0017 m

 

 

162

APPENDIX.B

Subprogram Unit IDDEL

This subprogram handles the calculation of the predicted values
using the user supplied set of differential equations. It also calcu-
lates the matrix of sensitivity coefficients. In support of this
routine the user must supply a subroutine called DEQIDD in the form:

SUBRWTINE DEAIDD(V.N8.TIB.WDT.INSUT.BETA)

REAL BETA(S).DNSDT}DVDT.NS.TIIE.V

OPTIONAL COIIANDS DBPENDING ON’NEEDH
. INTEGER IERUG.AOONFL
REAL CISTEP. CIZERO. CSSTEP. CSZERO. RADIUS. VIEAD

OOIImN/BUG/DEBUG

COINDN/BLOCIS/RADIUS

OOHMDN/BLOCK6/CSZERO.CSSTEP.CIZERO.CISTEP

OOHNDN/BLOCR9/VDBAD

OOIION/BLOCIOIACONFL
An example of the type of routine required is contained in Appendix D
of the present work. As illustrated in this example the user can sup-
ply additional support routines such as the one shown which calculates

the concentration of the sample region using Equation 4.1.16 of the

present work. The input variable are defined in these routines.

Note that this routine can be made to use a one-directional
derivitive (forward or backward difference) or a central difference
chosen via the variable CDIFF. In the present version the routine

always uses the central difference unless the parameter is at the

163
value of one of the user specified bounds. If this occurs then the
routine will use a one directional difference approximation in the
direction towCTds the center of the bounded region. The statements
needed to change this to a selectable option are included in the pro-
gram but in the current version are nulled out via a comment

indicator.

0001
0002

0003
0004
0005
0006
0007
0008
0009

0010
0011
0013
0014

0015
0016
0017
0018
0019
0020
0021

0022

OOOOOOOOOOOOOOOOOOOOGOO

COO

0060

000

164

THIS SUBRWTINE EVALUATES THE Vim-TIDE HISTORY FOR A GIVEN
SET OF PARALETER VALUES AND EXPERIDENTAL (ONDITIONS.

ITALSO

EVALUATES THE MATRIX 0F PARTIAL DERIVITIVES WITH RESPECT TO EACH
OF THE THREE PARAMETERS USING A CENTRAL DIFFERDICE APIROXIMATION

10

IF AND ONLY IF XFLAG=1 IN THE RWTINE CALL.

SUBTOUTINE INPUT IS:

N INTEGER NUBBER OF DATA POINTS
TIDE REAL(SO) VECTOR OF INDEPENDENT VARIABLE VALUES
BETA REAL(3) VECTOR 0F PARABETER VALUES:
DEFINED BY SUBRGJTINE DEQIDD
XFLAG INTEGER CONTRG. VARIABLE FOR SUPRESSION OF

DERIVATIVE CALCULATION .

SUBRGJTINE OUTPUT IS:

ETA REAL(SO) VECTOR OF CALCULATED IEPENDENT VARIABLE
VALUES
X REAL(50.3) MATRIX 0F PARTIAL DERIVITIVES
PRwRAIBD BY: STEVE NWLEV
4/ 83

SUBRGJTINE DDDEL( N. THE. ETA. X. BETA. IBBETA. IBBETA. XFLAG)
DIIENSION TIBEUO) . ETA(50) .X(50 .3) .BETA(3)
.BETAD(3) .DELB(3) .IBBETA(3) . IBBETA(3)
INTEGER XFLAG. N. DEBUG. CDIFF. P
REAL NS. NEERO. IBBETA
COMIDN/ BIB/ DEBUG
COMDDN/ BIOCKB/ DELB
COMDDN/ BLOCK4/ DT
COMIDN/ BLOCK6/CSZERO. CSSTEP. CIZERO. CISTEP. DPERMS. DPERMI. DEPTH
COMIDN/ BLOCKS/P. CDIFF

SET INITIAL mNDITIONS

NEEROBO .0

IF(CSZERO. NE.0 .0) NEERO-CSZERO/CSEXT (0 .0)
Y1-1 .0

NS-NEERO

CALCULATE PREDICTED VALUES FROM
USER mm. FOR GIVEN PARAIETER VALUES

CALL RK4(0.0.TII_E(1).Y1.NS.UT.BETA)
ETA(1)=Y1
NN-N-l
DO 10 I=1,NN
CALL RK4(TIE(I).TIB(I+1).Y1.NS.UT.BETA)
ETA(I+1)-Y1
CONTINUE

GECK IF DERIVITIVES ARE Nab

IF(IFLAC.NE.I)RETURN

0024

0025
0026
0027
0028

0029
0030
0031

0032
0033
0034
0036

0038
0040

0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054

000000

6000

15

COO

00000 000

GOOD

000

18

20

25

165

IF CDIFF.EQ.1 THEN CALCUTATE X
USING CENTRAL DIFFERENCE

IF CDIFF.NE.1 TEEN USE FORWARD
DIFFERENCE

[=2
IF(CDIFF.EQ.1)K=2

ZERO OUT DERIVITIVE MATRIX:

DO 15 131.3
DO 15 I81,N
X(I.J)-0.0
CONTINUE

SET UP DUMMY PARAMETER VECTOR

DO 18 I=1.P
BETAD(I)=BETA(I)
CONTTNUE

EVALUATE dV()/dBETA()

DO 30 LILP
DP‘0.01'ABS(DELB(L))
IF(DP.ED.0.0)DP-0.001‘BETA(L)
IF(L.EQ.3.AND.DP.EQ.0.0)DP*0.001

CHECK TO SEE IF CURRENT'PARAMETER
IS AT ONE OF PARAMETER BOUNDS AND
ROUTE ACCORDINGLY

IF(BETA(L) .LE.IBBETA(L) )GOTO 27
IF(BET‘A(L) .GE.IBBETA(L) )GOTO 26

THIS SECTION FOR.FORIARD OR CENTRAL
DIFFEREHE CALCULATION AS SET BY CDIFF

DO 25 181,!
BETAD(L)=BETA(L)+DP
Y1=1
NSBNSZERO ‘
CALL RX4(0.0.TTME(1).Y1.NS.DT.BETAD)
X(1.L)!X(1.L)+((XI-ETA(I))/DP)/X
DO 20 I-1.NN
CALL RK4(TIIB(I) .TIIE(I+1) .Y1.NS. UT. BETAD)
X(I+1.L)-X(I+1.L)+((Y1-ETA(I+1))/DP)/K
CONTINUE
DPB-DP
CONTINUE
GOTO 29

THIS SECTION FOR SPECIAL CASE
REQUIRING FORIARD 0R BACXIARD

0055

0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066

0067
0068
006 9
007 0

C
C
C

26

27

28

29
30

166
DIFF IF AT PARAIETER BWNDS

CONTINUE
IF AT UPPER LIMIT USE BACK-DIFF.
DP=-DP
CONTINUE
BETAD(L) =BETA( L) +DP
Y1=1 .0
NS=NEERO
CALL RK4(0.0 .TIDE(1).Y1.NS. D'T. BETAD)
X(1.L)=X(1.L)+((Y1-ETA(1) )lDP)/K
DO 28 I=1.NN
CALL RX4(TIBE(I) .TIDE(I+1) .Y1.NS. D'I'. BETAD)
x<1+1 . L) =X(I+1 . L) +(Y1-ETA(I+1) ) IDP
CONTINUE

HID SPECIAL SECTION

BETAD ( L) =BETA( L)

CONTINUE
RETURN
m D

167

APPENDIX C

Subroutine RI4

Subroutine RX4 is e nulericel interretion routine for two sinul-
tenious ordinery differentiel equetions. It is besed on e verietion
of the fourth order RunIe-Kutte nethod developed by Gill. [27]. end

presented in deteil by Ronenelli. [28].

The user inputs the velues of the lower end upper linits of
intezretion. T1 end T2 respectively. the velues of the two dependent
veriebles et the lower linit. Yl end Y2, end the epprorinete size of
the tine step increment to be used. The progren will celculete the
ectuel number of integretion steps fro-.11 to T2. NSTEPS. since the
user specified increlent ney not divide evenly into the totel tine
(T2-T1). The routine will then eveluete the velue of the dependent
veriebles et the tine T2 end plece these velues in 11 end Y2 before
returning to the ceiling routine. All other input veriebles ere

unchenged.

The routine requires en externel user supplied subroutine nmled

DEQMDD of the forn:

SUBRGITINE DEQDDUI .Y2 . TIE. DYlDT. DY2DT. BETA)

 

168

This subroutine should eveluete the derivitives of the dependent veri-
ebles. Y1 end Y2. with respect to the independent verieble. TIE,
besed on the current velues of Y1. Y2. TIME. end the vector of perene-
ter velues BETA. These velues ere then pieced in DYlDT end DY2DT es
the subroutine output. The velues of Y1. Y2. TIE. end BETA should

not be eltered in the routine.

Subroutine K4 conteins two debugging stetenents which will fleg
the cell to end successful return from the user supplied subroutine

DMDD. The level for output of these flegs is DEBUG.GE.1.

0001

0002
0003

0004
0005
0006

0007
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023

OGOOOOGOOOOOOOOOOOOOOOOOD

0060600

169

SUBRGITINE RK4(T1 .12 .Y1 .Y2 .DT. BETA)

TRIS SUBRGJTTNE IS A RUMB-KUT'TA NMRATION ROUTINE OF THE
FORTE ORDER FOR TWO SIELTANIOUS ORDINARY DIFFERWTIAL

'EE LIMITS OF INTEGRATION
TEE INITIAL FUNCTION VALUES AT T].
APPROXIMATE TIE INCREDENT (STEP SIZE)

EQUATIONS.
INPUTS TO THE RWTINE ARE:
T1.T2 REAL
'Y1.Y2 REAL
DT REAL
BE'TA

REAL(3) VECTOR OF PARAMETER VALUES

OUTHJT FROM THE RWTINE IS:
Y1 AND Y2 ARE REPLACED BY THE NEW VALUES AT TIME T2
T]. AND T2 ARE UNmANGED
DT IS UNmANGED

EXTERNAL FUNCTION mummm ARE FOR A SUBRWTDIE OF THE FORM:

SUBRCIITINE DEQIDD(Y1 . Y2 . TIE. DYlDT. DY2DT)

THIS EXTERNAL SUBRGITINE EVALUATES TEE ERIVITIVES 0F Y1 AND Y2
WITH RESPECT TO TIE AND RETURNS THEM IN DYlDT AND DY2DT

REPECTIVEY.

TEE VALUES 0F Y1. Y2. AND TIE SEWLD NOT BE

EDIFIED IN THE RWTINE.

DIWSION A(4) .B(4) .C(4)
REAL A. B. C, DT. DYlDT. DY2DT, B, QX. Q]. . m .X. Y1.Y2

1 .BETA(3)

INTEGER FLAG. I. J. NST'EPS. DEBUG

COMEN/ BIB/ D-m

DATA FLAG/ 0/

NOTE THAT FOR SPEED OF EXECUTION IN REPETED CALLS TO RX4 TEE
VALUES OF A. B. AND C ARE EVALUATED ONLY THE FIRST TIE TERWGE
TEE RWTINE. THEY WILL REMAIN UNCEANGED TRERAFTER.

IF FIRST TIE TERGIGR EVALUATE A. B.

ANDC

IF(FLAG. NE.0) GOTO 10

A(l)=0.5
SR-SQRT(0.5)
- A(2)=1-SR
A(3)=l+SR
A(4)=1.0/6.0
B(1)-2 .0
B(2)-1 .0
B(3)=1 .0
B(4)=-2 .0
C(1)=0.5
C(2)-A(2)
C(3)-A(3)
C(4)=0 .5
FLAG-l

10 CONTINUE

0024
0025
0026
0027
0028
0029
0030
0031
0032
0034
0035
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
. 0048

170

EVALUATE ACTUAL STEP SIZE

APPROX STEP SIZE. DT.

MAY NOT FIT IN'ERVAL

EVE‘JLY SO EST RECALCULATE

30
20

901
902

NSTEPSIIINTT (T2-T1) /DT)
E-(T2-Tl)/FLOAT(NSTEPS)
X==T1
DO 20 I=1.NSTEPS
QX=0.0
0180.0
(112-0.0
DO 30 J=1,4
IF(DEBUG.GE.1)TYPE 901
CALL DEQDDD (Y1 . Y2 . X. DYlDT. DY2DT. BETA)
IF(D.UG.GE.1)TYPE 902
X=X+E*A(J)‘(1.0-B(J)‘QX)
QX=QX+3.0‘A(J)‘(1.0-B(J’)‘QX)-C(J)
Y1=Y1+E‘A(J)‘(DY1DT-B(J’) '01)
01=QI+3.0'AU)‘(DYlDT-BUPQD-CUHDYlDT
Y2-Y2+E'A(J) ‘ (DY2DT-B( J') ‘02)
02:02+3 .O‘AU) ‘ (DY2DT-3(1) ‘QZ)-C (J’) ‘DY2DT
CONTINUE
CONTINUE
RETURN
FORMAT(' CALLING DEQIDD')
FORMAT( ' BACK FROM DEQIDD')
DID

171

APPENDIX D

Subprogrem Unit DEQMDD end CEXT

The routines described here ere exemples of the routines which
the user must supply when instelling e new trensport nodel for use
with MARBOX. MODEL. end 214. These exemples ere elso vieble for cel-
culetion purposes. The model is thet of ledem end Ketchelsky
formuleted for binery flow problems (see Chepter 3 Section 3 of the
present work). It is e very simple problem to nodify this model for
the cese of osmotic shrinkege only by setting ud-O end c-1.0 et the

top of the routine.

The support routine CEXT uses Equetion 4.1.16 to celculete the
concentretion of the semple region besed on the perneebility of the
dielysis membrene end the thickness of the region. CSEXT is for the
permeeble solute (reletive to the cell or liposome membrene) end CIEXT
is for the impermeeble solute (egein reletive to the cell or liposome

membrsne).

Note thet it is in routine DEQMDD thet the definitions of eech of

 

 

 

172

the peremeters in the vector BETA ere defined reletive to the modeling
equetions. As fer es MARBOX is concerned the peremeters ere not
defined implicitly. The user chooses which pereneter will be essigned
to BETA(l) end which to BETA(2) etcetere. Thus the user must neintein

consistency between the input end the modeling routine.

 

 

0001

0002

0003
0004
0005
0006
0007
0008

0009

0011
0012

0013
0014
0015

0016

0017

OfiGOOOGOOOOOOOOOOOOOOOOOOOOO

0006

173
SUBRWTINE DEQED(V. NS. TIE. DVD'T. INSDT. BETA)

THIS SUBRGITINE EVALUATES THE K-K SET OF PERMEABILITY
EQUATIONS IN DIFFERRTIAL FORM FOR THE VALUES 0F DV’IDT
AND INS‘IDT WHERE WE STAR IMH..IES NWDIDENSIONAL VALUES.

NOTE THAT IN THIS VRSION OF THE RWTINE mE TIE FACTOR
HAS NOT BER NWDIDENSIONALIZED. THIS LEAVES THE EQUATIONS
IN A PARTIALLY NWDITENSIONAL FORM ONLY.

THIS VERSION THEREFORE RFDUIRES THE INITIAL RADIUS OF
THE (ELL TO BE PASSED TERGIGH (OMEN BLOCK 5.

THIS RCIITINE ALSO WNTAINS THE FLAG FOR CONSTANT AREA
ASSUMPTION

INPUT TO THE RWTINE ARE THE VALUES OF N. AND V‘ AS WELL
AS THE TIE AT WHICH THE SEUTION IS DESIRED.

NOTE THAT ma TIE IS ONLY RRUIRED TO EVALUATE TEE
EXTERNAL SG.UTION (ONCENTRATION.

NOTE ALSO THAT THE K-K SET USED INCLUDES THE TERMS TO
HANILE IMPEREABLE SQUTES AS REL.

IT SHOULD ALSO BE NOTED THAT IT IS IN THIS RGJTINE THAT
m3 PARAE'ERS ARE DEFINED. IN THIS VRSION:
BETA(l) = Pw (80..va PERMEABILITY) (CM“4/MOLE-SEC)
BETA(Z) = Ps (S(LUTE PEREABILITY) (CM/SEC)
BETA(3) = SIGMA (REFLECTION (OEFFICIENT)

REAL BETA(B) .NS. Pl. PS. SIGMA

sum IS AN INTER DEBUGING VARIABLE
WITH ELT‘Y LEVH. CAPABILITY

INTEGR EBUG. ACONFL
COMlDN/ Bm/ DEBUG
COMIDN/ BLOCKS/ RADIUS
COMEN/ BLOCK6/ CSZRO. CSSTEP. CIZRO. CISTEP. DPERMS . DPERMI. DEPTH
COMIDN/ BLOCK9/ VDEAD
(1)le BLOC10/ ACONFL
CHECK FOR NFBITIVE VEUE
IF(V. LT.0 .0) GOTO 99
EV ALUATE EXTERNAL CONCENTRATIONS
CSWT-CSEXNTIE)
CIOUT‘BCIEH (TIE)
SET UP PARAETERS FROM BETA
P'=BETA( 1) ‘18 .015‘18—4
PSBBETAQ) '1E-4
S IG MABBETA( 3)
EVALUATE DV/DT
DVDTt-(PI‘S .0/ RADIUS) ‘
1 (SIGMA‘CSGIT‘(1.0-(NS/ (V-VEAD) ))
1 +CIOUT-(CIZRO‘(1-VDEAD) / (V-VIEAD) ) )
CHECK FOR Nm-CONSTINT AREA
IF (ACONFL.m.1)DVDT-DVUT‘(V”(2 ./3 .))

 

0019

0020
0022

0023
0025
0026

0027
0028
0029
0030
0031
0032

0033

99

901
902
903
904

2

EVALUATE IN S/ DT

174

X1=( (1 .0-(NS/ (V-VIEAD) ) )‘PS‘3 .0) /RADIUS
CHECK FOR NOV-CONSTANT AREA

IF (A(DNFL.HI.1) X1=X1‘ (V”(2 .l3 .))
DNSDT=X1+(DVDT‘(1 .0-SIGMA) " (1 .0+(NS/ (V-VIEAD) ) ) l2 .0)

IF(D-UG. GE.2)PRINT 903 . DVDT. INSDT. CSGIT

RETURN

CONTINUE
nus SECTION FOR NRATIVE VEUE
PRINT 904.“. PS. SIGMA

STOP

FORMAT( ' CALLING CSEKT')

FORMAT(' BACK FROM CSEH')

FORMAT(' DVUT =- '.m5.7.'
FORMAT(' ”‘ ERROR CONDITION IN DEQIDD ”"J.
1 ' NRATIVE VEUE RmUNTERED. PARAETER VALUES: '.I.

' PI

RD

3 'pmSer/p' PS

INSET = '.E15.7.'

. '3m5e63/D'

CSGJT = ' .F10.4)

SIGMA 8 '.E15.6)

 

0001

0002
0003
0005
0006
0007
0008
0009
0010

0001
0002
0003
0005
0006
0007
0008
0009
0010

000600000000

1

C

0

175
FUNCTION CSEXNTIE)

THIS FUNCTION IS A USR GRRATED FUNCTION TO EVALUATE HE
EHERNAL S(LUTE (DNENTRATION AS A FUNCTION OF TIE.

THIS IS OF PARTICIILAR INTREST IN THE DIFFUSSION (EAER WHRE THE
EXTERNAL SQUTION CANNOT BE ASSUMED TO UNDERGO A STEP CHANGE.

NOTE THAT IF THE INVESTIGATOR BEIEVES THE STEP (RANGE IS
APPROPRIATE EARLY HAVE THIS FUNCTION RETURN A CONSTANT VALUE.
IN THIS VRSION THIS IS NNE IF mE DIALYSIS EERANE
PEREABILITY IS GREATR THAN 998.0 (CM/SEC).

COMEN/ BLOCK6/ CSZRO. CSSTEP. CIZRO. CISTEP. DPERMS. DPERMI . DEPTH
IF(DPERE . LT.998 .0) GOTO 10
CSEXT=CSSTEP
RETURN
CONTINUE
CSEXT=CSSTEP+ ( CSZRO-CSSTEP) I'EXP ("DPERMS 'TIE/ DEPTH)
RETUR
RD

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

10

FUNCTION CIEKT (TIE)
COMle BLOCK6/ CSZRO. CSSTEP. CIZRO. CISTEP. DPERIB. DPERMI. DEPTH
IF(DPERMI . LT.998 .0) GOTO 10
CIEXT=CISTEP
RETURN
CONTINUE
CIEXT¥ICISTEP+ ( CIZRO-CISTEP) l'EXP (-DPERMI ‘TIE/ DEPTH)
RETURN
RD

 

 

 

APPENDIX E

Subproyrm Units TAKEP'I‘ and CIRCLE

 

 

 

 

176

0001

0002

0003

0004
0005
0006
0007
0008
0009
0010
0011

OOOOOOOOOOOOOOOOOOOGOOOGOOOOOOOGOOOOOO

PRCBRAIBD BY:

mun. MD1Y

ENDZX. WDZY

II

X(1) . Y(1)

1(2) .Y(2)

X(II) . Y( II)
WHERE:

E‘IDlX INTEGER

MDlY INTEGER

mm INTEGER

ENDZY INTEGER

II INTEGER

X INTEGER

Y INTEGER

171
THIS IRMRAH UNIT INTEGRATES A SERIES OF ARCS DEFINED
BY THREE POINTS EACH. THE FINAL AREA IS THE SUM OF THE
INTEGRATED PARTS DIVIDED BY TWO. TRIS PROCEWRE IS
DOGHENTED IN CHAPTER 4 SECTION 4 OF THIS THESIS WORK.

INPUT TO HIS RWTINE IS THRIIJGE A FILE NABD POINTS.DAT
THIS FILE SEWLD CONTAIN IN THE FOLLWTNG ORDER:

‘ X COORDINATE OF FIRST SCALE WDNNT

Y COORDINATE OF FIRST SCALE DIDPOINT
X COORDINATE OF SEmND SCALE EJDPOINT
Y COORDINATE OF SEWND SCALE WDNNT
NUDBER OF DATA POINTS

SET OF I COORDINATES

SET OF Y COORDINATES

TEE PRmRAM OUTPUT IS WRITTEQ INTO A DATA FILE NABD
POINTSJDT TEE PRCBRAN WILL DESTORY ANY EXISTING FILE
WITH THIS NAB DURING EXECUTION SO TEAT ANY FILE TO BE
SAVED MUST BE REWARD BEFORE RUNNING CIRCLE AGAIN.

STEVE NWLEW
3/ 83

DIWSION 1(20) .Y(20) .XX(3) .YY(3) . SUBA(3.3)
1 .SUBB(3.3).SUBC(3.3).COEF(3.3).R(20).THETA(3)
INTEGER LJ’, II,P,K(RINT,L. SFLIGJI

_, D-UG

COMIDN/ BUG/ DEBT];

COOOOOOOOOOOOOOOIQOOOQOOOOOOO

C

C DATA INPUT SECTION:

C

OPEN(UNIT‘-2.NAlE-'DK:POTNTS.DAT' .TYPEB'mD')
READ (2.‘)IDJD1X.IEND1Y

READ (2.’)IE~W2X.IE€D2Y
F1-ABS(FLOAT(IE{D1X-IENDZX) )
F2-ABS(FLOAT(IEND1Y-IEWD2Y))
SCALE=SIRT(F1"2+F2”2)

READ (2.‘)N
DO 10 I-1,N

 

 

0012
0013
0014
0015
0016
0017

0018
0019
0020
0021
0022
0023
0025
0027
0029
0031
0032
0033
0034

0035
0036
0037
0038

003 9
0040
0041
0042

0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060

10

0000000

’19

29

22

23

178

READ (2.‘)IX.IY

X(I)=FLOAT(IX)

Y(I)==FLOAT(IY)
CONTINUE
CLOSE(UNIT‘-2 .DISPOSE= ' SAVE')
OPEN(UNIT‘-4.NADE='DK:POINTS.WT' .TYPEB'NEW')

BHiIN CALCULATION PROCEDURE:

omnttmutt TOP 012 MAIN mmm LOOP an.

XMAX=0.0
XMIN=IOO0.0
YHAX=0.0
YMIN=IOO0.0
DO 19 I=1,N
IF(X(I) .G'T.XHAX)XMAX=X(I)
IF(X‘I).LT.XMTN)XMIN=X(I)
IF(Y(I).GT.YMAX)YMAX=Y(I)
IF(Y(I).LT.YMIN)YNIN=Y(I)
CONTINUE
xo-(XMAX+XMIN)/2.0
YO=(YMAX+YUN) l2 .0
DO 20 I=1,N
CONVERT TO RADIAL
R(I)=SQEIT(X(I)-XO)“2.+(Y(I)‘YO)“2.)
CONTENUE °
ARE!F0.0
DO 21 1'1,N
ROTATE COORDS AND CALC TEETA
TEETA‘I)‘ATAN((Y(I)-YO)/(X(I)’XO))
TEETA(2) =ATAN( (Y(I+1)-YO) ’(X(I+1)"XO) )‘THETAUJ
TEETA(2)=ATAN((Y(I+2)-YO)/(X(I+2)-XO))‘THETA(1)
TEETA(1) =0 .0
CALCULATE A. B, AND C:
DO 22 131.3
COEF(J.1)'1.0
COEF(J.2)'TEETA(J)
COEF(J.3)=EEETA(J)“2.
CONTINUE
DETXF(COEF(2.2)‘COEF(3.3))‘(COEF(2.3)‘COEF(3.2))
DO 23 151.3
SUBA(J:1) =R(I+J"1)
SUBA(L2)'wEF(J.2)
SUBA(L3)=COEF(L3)
SUBB(J.2)=E(I+J-1)
SUBB(J,3)=COEF(J.3)
SUBC(J.1)'(X)EF(L1)
SUBC(J.2)'COEF(J.2)
SUBC(J:3).R(I+J’1)
CONTINUE
AFDETERI“ SUBA) IDEn

0061
0062
0063
0064

0065

0066
0067
0068
0069
0070
0071
0072
0073
0074
0075

21

C

905

906

907

1

179

B=DETERM(SUBB)/DETX
C=DETERM(SUBC)/DETX
ThTEETA(3)
AREA.- AREA + 0.5"I ( T‘(A“2.) + A!B‘(T"2.)
+ (2.‘A‘C+(B“2.))‘(T‘*3.)/3.
+ B'C*(Tu4.)/2. + (C"2.)‘(T"5.)/5.)
CONTINUE
END PROGRAM SEQUENCE
WRITE(4.905)AREA
FORMAT(' AREA:- '.E15.8.'
PI=4.0‘ATAN(1.0)
RADIUS=SQRTTAREAIPI)
WRITE(4.906)RADIUS
FORMAT(' RADIUS = '.F10.4.’ UNITS')
WRITE(4.907)SCALE
FORMAT(' SCALE I 1 :
STOP
END

UNITS')

'.F10.4)

 

 

 

0001
0002
0003

0004
0005

1
2

180

FUNCTION DETERM(X)

DILENSION X(3.3)
DETERM-X(1.1)’(X(2.2)‘X(3.3)-X(2.3)"X(3.2))
- X(1.2).(X(2.1)‘X(3.3)-X(2.3)‘X(3.1))

+ 1(1.3)‘(X(2.1)'X(3.2)-X(2.2)*X(3.1))
RETURN
END

 

 

0001

0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022

OOOOOOOOOOOGOOOOOOOOOOOGO

181
WWW POINTS

THIS PRMRAM UNIT IS LOCAL TO ma IMAGE ANALYSIS
LABORATORY OF MIGIGAN STATE UNIVERSITY. ITS PURPOSE
IS TO READ IN A SERIES OF X-Y COORDINATE POINTS. THE
USER SIGNIFIES THAT THE (IIRSOR IS ROPERLY POSITIONED
FOR THE NEXT POINT BY DEPRESSING TEE BACK-SPACE KEY ON
THE (ONS(LE TERMINAL. TEE PRCBRAM WILL mm RING TEE
TERMINAL BEL WHEN READY FOR ANOTHER POINT.

nus PRmRAM HIST BE LINKED TO THE LIBRARY TVLIBJBJ'
ON THE SYSTEM VGJIIB. THIS LIBRARY CONTATINS TEE EXTERNAL
ROUTINE CURSOR. TEE PRmRAM IS CURRENTLY SET UP TO READ
FIRST:

TWO SCALE EWDPOINTS FOR REFERENCE SCALE

20 DATA POINTS ARGIND mE ENE OF A CIRCULAR IMAGE

OUTPUT FROM THE RWTINE consrsrs OF THE DATA POINTS IN
IN'TFfiER FORMAT (IS) AND HE NUDBER OF DATA POINTS (20) IN
PROPER FORMAT FOR DIRECT EJTRY INTO THE RWTINE CIRU..E.FOR
TEE OUTPUT FILE NAE IS POINTS.DAT WEIG IS THE INHIT
FILE NATE FOR RWTINE CIRCLE.FOR.

PRCBRAIED BY: STEVE NWLEQ
319

INTWER I.K.X.Y
EHERNAL OWLCURSOR
OPHT(UNIT‘-3,NA|E='DK:POINTS.DAT'.TYPEP'NH')
TYPE 905
FORMAT(' mm SCALE WU POINTS NW')
CALL WRSOR(X.Y)
WRITE(3.901)X.Y
FORMAT(ZIS)
CALL WRSOR(X.Y)
WRITE(3.901)X.Y
'RITE(3.902)
FORMAT( ' 20 ')
TYPE 903
FORMAT(' ENTER 20 DATA mIN'TS NW')
DO 10 1.1320
CALL CURSOR(X.Y)
WRITE(3.901)X.Y
CONTINUE
CLOSE(UNIT‘=3 .DISPOSE= ' SAVE')
STOP
m0

905

901

902

903

10

 

 

APPENDIX F

Subprogram Unit CURGEN

This program unit is for use in generating predicted curves from
specific parameter values independent of the parameter estimation
routine. The routine requires the external routines RK4 and DEQMDD as
described in the previous sections. The input-output format is des-

cribed in the program listing.

182

0001
0002
0003
0004

GOODOOOOOGOO(50000606000OOOOOOOGOOOOOOOOOOOOOOOOOOOO

183

THIS ROUTINE IS A DRIVER FOR GENERATING SEUTIONS
TO A SET OF DIFFERENTIAL EQUATIONS FOR THE mANSPORT
OF MATERIALS ACROSS A ELL EBBRANE. THE RETINE
SEVES TEE IDDELING EQUATIONS FOR A TIE SEUTION OF
VEUIE AND SEUTE ENTEWT.

INPUT TO THE RETINE IS THREGE LEICAL UNIT 3 SO
EAT UNER RT-ll TEE INPUT FILE SHELD BE NADBD

"FTN3.DAT".
INPUT IS IN THE FREE FORMAT FORM AND IN THE ORDER:

D-UG

DTINT. DTOUT. TMAX
CIZERO. CISTEP
CSZERO. CSSTEP
DPERMI . DPERBB. DEPTH
PI , PS, SIGMA

RADIUS

VEAD

AENFL

E

D-m INTEER VARIABLE LEVE. -113 ETPUT
(DNTRG. VARIABLE

DTINT REAL STEP SIZE FOR INTEGRATION (SEC)

DTET REAL INTERVAL FOR PRINTING OF
VEUIE AND SEUTE ENTWT (SEC)

TMAX REAL LEIGTE OF TIME TO INTERATE OVER (SEC)

CIZERO REAL INITIAL IMPERMEABLE SALUTE
ENENTRATION

CISTEP REAL NH IMPERMEABLE SEUTE CONC.

CSZERO REAL INITIAL EREABLE SEUT'E ENC.

CSSTEP REAL NEW ERMEABLE SEUTE CONC.

cease NOTE ALL ENC. IN UNITS HOLES/CC) ””‘
DPERMI REAL DIALYSIS EBRANE PEREABILI'TY
TO IMPERMEABLE SEUTE (CM/SEC)

DPERMS REAL DIALYSIS MEDBRANE EREABILITY
TO PERMEABLE SEUTE (CM/SEC)

DEPTE REAL DEPTH OF SMPLE REIGN IN DIFFUSION
EAIBER (CM)

P" REAL HYDRAULIC OR 8(1va (WATER)
PERMEABEITY OF CELL (MIECNS/ SEC)

PS REAL SEUTE PERMEABILITY OF CELL (MIEWS/ SEC)

SIGMA REAL INTERACTION TERM

RADIUS REAL INITIAL ELL RADIUS (CM)

VEAD REAL FRACTION OF INITIAL VEUE EVOTED

OSDDTICALLY INACTIVE VEUE
A(DNFL INTEGER FLAGS FOR AREA AS FUNCTION OF VEUIE
ACONFL I 0 IMHJES ENSTANT AREA
AENFL - 1 IMHJES AREA-FUD

DIIENSION BETA(3) .VSTORE(50) ,TIEUO)
REAL NS. Pl. PS. SIGMA

INTEGR sum. ACONFL

COMDN/ BLOCKS] RADIUS

 

0005
0006
0007
0008

0009

0010
0011
0012
0013
0014
0015
0016
0017
0018

0019
0020
0021
0022

0023
0024
0026
0028
0029
0030
0031
0032
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0045
0046
0047
0048
0050
0051
0052
0053
0055
0057
0058
0059
0060
0061
0062

10
15

184

COMIDN/ Ems/csznno, CSSTEP, CIZERO. CISTEP. DPERMS, DPERMI. DEPTH
COMWN/BwCK9/VDEAD
comm BLOC10/ ACONFL
COMIDNIBUG/D-UG
OPEN FILE FOR OUTPUT
omN (UNIT=6,NADE='DY1:OUTPUT.DAT".TYPEs'NEW')
DATA ENTRY BLOCK
READ(3."')DEBUG
READ(3.‘)DTINT,UI‘GJT.TMAX
READ(3 . ‘)CIZERO. CISTEP
mn(3.t)cszmo,cssmp
READ(3.‘)DPERMI.DPERMS.DEPTH
READ(3.')PW,PS,SIGMA
READ(3.‘)RADIUS
READ(3.‘)VDEAD
READ(3.")A(I)NFL
ECHO OUTPUT OF EXPERIMENTAL UONCITIONs
mun 6 . 905) CIZERO. CISTEP. cszmo. CSSTEP
“mu 6 , 906) DPERMI , OPEENS. DEPTH
mTE<6.904)Pw. PS. SIGMA
WRITE(6.907)RADIUS.VDEAD
BEGIN CALCULATION PROCEDURE
NSTEPS=mT(TMAX/DTOUT)
IF(NSTEP.CT.50)wnrrE(6.911)
IF(NS‘I‘EPS.GT.50)NSTEPS=50
BETA(1)=PI
BETA(2)-PS
BETA(3)-SIGIIA
Ns=0.0
IF(Cszm0.NE.O.0)Ns-cszm0/CSEIT(0)
TIIE2=0.0
v-1.0
WRITE(6.900)
WRITE(6.901)TIIE2.V.NS
IFLAG-O
ROUND-0
no 10 I=1,NSTEPS
TInEI=TnE2
TIDEZ=TIBB1+DTOUT
IF(DEBUG.GE.1)TYPE 903
CALL RK4(TIIE1 , TIIEZ . v. Ns. UTINT, BETA)
VSTORE(I)=V
TInE(I)=TIm2
IF(D-UG.GE.1)TYPE 902
IRITE(6.901)TI)E2.V.NS
CHANGE=ABS(VSTORE(I )-VS'mRE(I-1) )
KOUNTa-KOUNTu
IF(CBANGE. LT.0 .0001) IFLAG- mama
IF(IFLAG.GE.3)GOTO 15
CONTINUE
CONTINUE
WRITE(6.908) ‘
Vim-21.0
00 16 I-1.KOUNI'
IF(vsmnEm .LT.vm)vuIN=vsmEE(I)

 

 

0064
0065
0066
0067
0068
0069
0070
0071
0072
0073

0074
007 5
007 6

0077

007 8

007 9
0080
0081

0082
0083

16

20

901
900

902
903
904

905

906

907

908

909
911

185

CONTINUE

DO 20 I=1.KOUNT
VSTAB=(vstE(I)-VNIN)/(I.0-VMIN)
RADRAT=(ABS(VS'IORE(I)))”(1./3.)
wBITE(6.909)TIIE(I).VSTAB,BADEAT

CONTINUE

CLOSE(UNIT=6)

STOP

FORMAT(E15.4.IOX.£15.6.10X.E15.8)

FORMAT(/.' NOTE: NSTAR IS DEFINED AS IIIE NUDBER 0F MOLES'

1 ./.' PERMEABLE smUTE INSIDE DIVIDED BY THE EIT'EBNAL'

2 .l.’ PERMEABLE SILUTE mNCEN'IRATION AT TIDE t DIVIDED'

3 ./.' BY THE INITIAL vaNE'J

4 ./.sx,' TIDE (SEC)',9x.'va.UIE(t)/va.UME(0)',IOX.'NSTAB')

FORMAT(' BACK FROM n40

FORMAT(' CALLING RK4')

FORMAT‘(/.' HYDRAULIC PERMEABILITY - '.E15.6.' nmmS/SEC'

1 ,/,' SGJTI‘E PERMEABILITY = '.El5.6,' mmmS/SEC'

2 .I.' SIGMA . =- '.F6.4)

FORMAT(' INITIAL IMPERMEABIE 80.011! CONCENTRATION a '.E15.6.

I ' IDLES/CC',
2 .l.' NE! IIIBEBICIABLE SILU'IE CONCENTRATION - '.E15.6
3 ' IDLES/CC',
4 Jim INITIAL mBIIEABLE SOLUTE CONCENTRATION - '.E15.6
5 ' IDLES/CC',

6 .I.’ NE! BEBNEABLE 8mm CONCENTRATION = '.El$.6
7 ' mLES/CC')

FORMAT(/.' DIALYSIS IEIBBANE PEREABLILITYUJ.

1 ' TO INPEENEABLE SGJJTE - '.El5.6.' (CM/SEC)’.I.
2 ' TO PERMEABLE sum as '.E15.6.' (CM/SEC)’,/,

3 I.’ DEPTH 0F SAMPLE RHEION 8 '.E15.6.' (CM)')

FORMAT(I.’ INITIAL (ELL RADIUS = '.E15.6.' (CM)',/,

1 ' OSBDTICALLY INACTIVE V(LULE FRACTION 8 '.F6.4)

FORMATU ./ .SX. ' TIDE (SEC) ' .7X. ' (V(t)-Vmin)/(V( 0) -Vnin) ' .81.

1 'DSTAR=D(t)/D(0)')

FORMAT(ZX.m5.6.8X.E15.6.8X.E15.6)

FORMAT(/.' EXCEEDED MATRIX DIWSIONS. WILL MAKE 50 STEPS'J.)
DID

1.)

2.)

3.)

5.)

186

APPENDIX.G

Tabulated Osmotic Shrinkage Data

Experimental Conditions

All cells initially in 0.3 (camel/kg) saline solution.

For cells 1 and 2 the saline concentration increased
to 0.5 (osmol/kg). For cells 3 and 4 the saline con-
centration increased to 0.8 (camel/kg).

Initial radius as follows:
cell 1: r - 39.5 (microns)

cell 2: r - 40.4
cell 3: r - 38.0
cell 4: r t 38.5

Osmotic inactive volume as fraction of initial volume:
cell 1: 0.219
cell 2: 0.156
cell 3: 0.292
cell 4: 0.260

All experiments conducted at room.temperature using the
diffusion chamber. Dcpth of sample region - 228.3 (microns).
Cuprophan 80pm.msmbrane used.

 

Time
(sec)
0.0
10.0
20.0
30.0
40.0
50.0
60.0
70.0
80.0
90.0
100.0
110.0
120.0
130.0
140.0
150.0
160.0
170.0
180.0
190.0

200.0

cell 1
1.0

1.0

0.9479
0.8747
0.8651
0.8147
0.8417
0.7933
0.7752
0.7227
0.7170
0.7342
0.6836
0.6726
0.6863
0.6836
0.6863
0.6781
0.7001
0.7001

0.6974

187

V‘ = V/ v.

cell 2
1.0

0.964
0.903
0.847
0.791
0.729
0.752
0.758
0.730
0.706
0.693

0 .677

0.686 '

0.662
0.659
0.662
0.653
0.658
0.672
0.667

0.667

cell 3
1.0

0.9787
0.8320
0.7086
0.6416
0.6352
0.5961
0.5950
0.5604
0.5546
0.5575
0.5585
0.5575
0.5575

0.5575

cell 4
1.0

0.9703
0.8710
0.8096
0.7290
0.6250
0.5514
0.5595
0.5574
0.5615
0.5514
0.5374
0.5374
0.5314

0.5374

1.)

2.)

3.)

4.)

5.)

APPENDIX E

Thbulated Binary Flow Data

Experimental Conditions

All cells initially in 0.3 (osmol/kg) saline solution. New
solution also 0.3 (osmol/kg) saline content. .

New solution contains solute glycerol. For cell 5 the
concentration is 0.2 (osmol/kg). For all others the
concentration is 0.25 (camel/kg).

Initial cell radius as follows:
cell 5: r a 39.3 (microns)
cell 6: r 8 37.0
cell 7: r . 33.2
cell 8: r B 33.5

All experiments run at roan tamperature using the diffusion
chamber. CuprOphan type 80pm membrane used.

Depth of sample region in diffusion chamber a 228.6

(microns).

6.)

Osmotic inactive volume unknown in all cases.

188

 

 

Time
(sec)
0.0
10.0
20.0
30.0
40.0
50.0
60.0
70.0
80.0
90.0
100.0
110.0
120.0
140.0
160.0
180.0

190.0

cell 5
1.0

0.9765
0.9224
0.8196
0.7998
0.7255
0.6876
0.6764

0.6823

 

0.6571

 

0.6865
0.6764

0.7092

 

0.7192

w-

189

cell 6
1.0

0.9504
0.8106
0.7354
0.7064
0.6301
0.5991
0.5834
0.5941
0.5976
0.6178
0.6025
0.6102
0.6267
0.6315

0.6575

w vo

cell 7
1.0
0.98
0.96
0.94
0.90
0.848
0.83
0.779
0.729
0.713

0.713

0.681
0.697
0.666

0.681

cell 8
1.0
1.0
0.852
0.800
0.735
0.735
0.719
0.719
0.704
0.735

0.689

0.704
0.712
0.689

0.704

Time

(sec)
210.0
220.0

240.0

250.0'

270.0
280.0
300.0
310.0
330.0
340.0
360.0
370.0
390.0
400.0
420.0
480.0
540.0

600.0

cell 5

0.7498

0.7865

 

0.7861

 

0.8001

 

0.7873

 

0.8360

 

0.8686

190

V" = Vlvo

cell 6

0.6453

 

0.6999

 

0.7056

 

0.7379

 

 

0.7755

cell 7

0.666

0.697

0.666

0.713

0.697

0.729

0.713

0.745

cell 8

0 .719

0.719

0.735

0.735

0.766

0.800

0.784
0.835
0.870

0.888

 

BIBLIOGRAPHY

 

10

11

12

13

14

BIBLIOGRAPHY

Jacobs. M.B.. in Barron. E.S.G.. 112 W of, 11mm
“hm. University Press. Cambridge Msss.. 1952.

Zuethen. E.. and Prescott. D.M.. £913. M191. 3951151..
28(1953) .77.

Kedem. 0.. and Katchalsky. A.. 3193213133, E}, M119; 4951.
"Thermodynamic Analysis of the Permeability of Biological
Membranes to Non-Electrolytes". 27(1958) .229.

Onsagor. M 391.. 37(1931) .405; 38(1931) .2265.

Levin. R.L.. Mung BigLOgy. "Water Permeability of Yeast
Cells st Sub-Zero Temperatures". 46(1979).91-124.

Staverman. A.J’.. ML, T331, M 221223.91. Belg“ 70.344.

Kedem. 0.. and Ratchslsky. A.. .11! 1mm 21', 9:31:31 M.
"A Physical Interpretation of the Phenmenological .
Coefficients of Membrane Permeability". 45(1961).143.

Johnson. LA. . and lil son. T.A. . Logan}, .91 W21 Mon.

"Osmotic Volume Changes Induced by a Permeable Solute".
17(1967) .304-311.

Katchalsky. A.. and Curran. P.F.. W W
13 W: University Press. Cambridge Mass. 1965.

Stusnick- En Mi lull—0Q 13 2212111212: Mane.

Will-1! W. Ph.D Dissertation. State Univ.
of New York at Buffalo. (unpublished) 1971.

Stusnick. E.. and Burst. R.P.. M of, mm; mm.
"Numerical Determination of Membrane Permeability
Parameters". 37(1972) .261-271.

Pnpanok. “M!” In MI W 9.1 has; 8129192159. is

I» W Ens: £2.59 .LQ fl. Ph.D Dissertation.
Mass. Inst. Tech.. (unpublished) 1971.

Fenstermacher. J’.. and Johnson. J.A. . 31122.12 19.91381. 21 m.

211(1966).341.
Willbrandt. I” £11m Auk. 1.21.. m 254(1941).1.

1903.

 

15

16

17

18

19

20

21

22

23

24

26

27

28

29

30

191

Book. J.V.. and Blackwell. B.. glgtion M91; f9; Inger” Bgat
W M. (unpublished as of Aug. 1983).

Beck. J.V.. and Arnold. K.J.. Pgrggter Estimgtiog _iA Engineering
114 m. Wiley. New York. 1977.

Bard. Y.. M 283.9912; Estimation. Academic Press Inc..
New York. 1974.

Stewart. 0.11.. and Jacobs. 11.11.. M}, C911, 922. M2]...
1(1932).83.

Ligon. R.. Personal Communication. 1982.

Shabana. M.. Personal Communication. 1983.

Terwilliger. T.C.. and Salmon. A.K.. ml 21 £221.11 m.
77(1981) .549-570.

Lovelock. J.E.."Physica1 Instability of Human Red Blood Cells".
3.19.0112. 1... 60:692-696.

Mazur. P.. "The Role of Intracellular Freezing in the Death of
Cell Cooled at Supra-Optimal Rates". CM14(1977):251.

Rayliegh. J.'.S.. 111 1111221 21 m. Dover Publ.. New York.
. 1877 (reprinted 1945).

Dupuis. C.A... "A Preliminary Study of the Diffusion Through
Cuprophan Membranes of Sodium Chloride and Sucrose".
(unpublished) 1983.

Melkerson. M.N.. "The Determination of Dialyzing Membranes Effects
on Liposaue Response Characteristics". (unpublished) 1983.

Gill. 8.. Proc, We Philos, Soc. 47(1951).96-108.
Rcmanelli. M.J.. in Ralston. A.. M 312941.! 19; M

W. "Runge-Kutta Methods for the Solution of
Ordinary Differential Equations". Wiley. New York. 1960.

Tu. Sui-Ming. Personal communication. 8/83.

Jackowski. 8.. Leibo. S.P.. Mazur. P.. 119, m1 91 W1
1991.331. 212(1980):329-341.