BULLETIN OF MATHEMATICAL BIOLOGY

VOLUME 37, 1975

SIMULATION OF CELL B E H A V I O R : NORMAL AND ABNORMAL GROWTH

9 E. J. DAVISON Department of Electrical Engineering, University of Toronto, Canada

A mathematical model simulating a cell growing in a culture medium is obtained. Using this model, various behavioral patterns of the cell are obtained under different types of disturbances, in particular when (i) a Mg2§ deficiency experiment and, (ii) a split-dose ionizing radiation experiment are carried out, (iii) when disturbances on the rate constants of the biochemical reactions taking place in the nucleus of the cell are applied, and (iv) when the cell's interior components axe perturbed. The cell model results obtained agree well with experimental results for the Mg2 + and split dose experiments, and explain the mechanism of the split dose radiation experiment without the need to introduce additional axioms (e.g. healing processes) into the dynamics of the cell. Conditions are obtained which cause the cell to behave in a rapidly growing 'tumor-like' mode; it is shown that once the cell moves into this 'tumor-like' mode, its behavior is irreversible, i.e. if a disturbance of opposite type is then applied to the 'tumor' cell, the cell will not revert back to its original normal behavior.

1. INTRODUCTION The purpose of this paper is to find the possible mechanism which explains w h y a n o r m a l cell suddenly becomes 'tumor-like'. This will be done b y obtaining a m a t h e m a t i c a l model of a growing cell living in a culture medium, and then systematically examining the behavior of the cell when various disturbances are applied to different parts o f the m a t h e m a t i c a l model o f the cell. The m o t i v a t i o n for applying such a model approach is t h a t m a t h e m a t i c a l system techniques have been highly successful in explaining and predicting the behavior of technical-industrial processes (see e.g. Himmelblau and Bischoff, 427

428

E.J.

DAVISON

1968), and recently have given insight into the mechanism of various biological processes (see e.g. Mesarovic, 1968; Diamant et al., 1970; Palmby et al., 1974}. There has been to the author's knowledge no previous work done on this problem. The idea however of describing a cell using a mathematical model is not new. See for example Pollard (1960), Yeisley and Pollard (1964) who described a cell using 5 and 7 differential equations respectively. Heinmets (1964a), (1964b), {1966) has modelled a cell using 19 differential equations to study various transient responses of the cell assuming no cell division takes place, and his qualitative model shall form the basis of the model to be described in this paper. Glueck (1969) used Heinmet's model to extend some of Heinmet's results. Weinberg and Zeigler (1969) studied the transient behavior of a non-dividing cell using automata theoretic analysis. Other papers have been concerned with the modelling of tumor cells per Be, e.g. Jensson {1968) finds a model of tumor cells assuming that the cell population is tumorous, instead of finding a model of normal cells and then finding conditions under which these normal cells become tumorous. All the previous papers use either analogue or analogue-digital computers to study the mathematical model of the cell; this study shall use solely digital computers. In the synthesis of any model of a physical (or biological} system, certain assumptions as to what effects to include or exclude must always be made, and this is certainly true in the case of obtaining a mathematical model of a cell. It is also obvious that it is not necessarily clear (initially at least) what effects are more important and should be included, and what effects need not be included in such a model. In the ease of finding conditions for a cell to become 'tumor-like' however, one receives insight into its possible mechanism even if the model has not included enough effects; i.e. if the assumed model does not have the property of becoming 'tumor-like', then one can assume the mechanism of tumor-growth lies in those effects which have not been incorporated into the model. If the assumed model does have the property of becoming 'tumor-like', then one can assume the mechanism of tumor-growth lies in those effects which have been incorporated into the model. In the model to be described, any dynamics occurring outside the nucleus of the cell have been ignored, e.g. interaction effects of cells with neighboring cells have not been incorporated. It will be shown however that the present model is rich enough in structure to give 'tumor-like' behavior when certain disturbances occur. The paper is organized as follows: Section 2 describes the normal cell model, Section 3 studies the behavior of the cell model under Mg 9'+ deficiency, Section 4 deals with the case when the cell is subject to a split dose of ionizing radiation, and Section 5 studies those disturbances which cause the cell to become 'tumor-like'.

SIMULATION OF CELL BEHAVIOR

429

2. NORMAL CELL MODEL

This section deals with the qualitative behavior assumed for the model of a cell growing in a culture medium. 2.1 Normal Cell Chemistry. The topic of cell chemistry is a vast subject; some representative texts are given in Mitchison (1965), Teir and Rytomaa (1967), Mallette et al. ( 1971), and review articles in Medoff and Swartz (1967 ), Warner and Soeirs (1967). The basic components (states) of the cell system are: genes, messengers, templates, enzymes, repressors and activators. The following biochemical equations (Table I), taken from Heinmets (1966), are assumed to describe the cell's behavior. (See end of paper for definition of symbols used.) It is observed that all dynamic behavior outside the nucleus of the cell has been TABLE I

Complete Set of Reactions Occurring in the Cell 1.

k1 Ep + Pn ~ [EpPn]

14. B kl--~ 4 P. i 15, C kl-~P. i

2.

B' + E k-~2 B k. 3. GB + Pn ~ GB + B' k4 4. GC + Pn "~ Gp + .'I P k5 5. Gp + Pn -'g Gp + I'll) 6.

k6 GE + [EpPn] ~ GE + Ep + H

7.

k7 B + H --+ N P

9,

16. E k1-~6p. I. 17, E k~17 P. p I 18, P + E k~18 P. + E e

1

19. P. + E kI-~9 P + E i n 20. P. + E k-~20 P + E i a

P

21. Ep

C + Pa];-~ 9 [CPa]

I0. N + [CPa] ~

II. Np + [CPa] ~

ck~-~ k31 [EpC]

k22 22. E + Pi ~32 [Epi]

B§ C§ H§ E

B+Hp+C+E

+

P

23. Ep +

B k2~ k3 ~ [EpB]

12. Mk-~ 12 P n

13. Hp ~

Pn

24. N - ~ p .

1

25. N - ~ P . p I

26. p3. k'-~X

430

E. J. DAVISON

ignored in this model and that the extracellular nutrient pool Pe is assumed to be a constant (but its value can be varied). It is in addition assumed that the temperature of the cell is held fixed. Figure 1 gives a block diagram for the interconnections of the various components of the cell.

flNA iynthesis

ks f / ~ ~ / . ~ .~",1, /.;- . ~ I

I "~:Y,,'~

,).

J

km

_

] [

I

I

[

~/

/ 9 9

/

~

.

k=y I / I

~'~.~ t

, 4] /" 1 , 1 -

/ i l.~C- - ~

|

//I

I

/~I

|k=3 I !

T

INTERMEOIATE COMPLEXES

[r

II11 // I I

.......

IMzs**n=, ~oduci'~n

/lf

_@/

Pool Pi

I

',, / ~zl,i,!_/.

/

R*~nerafion

I:;' .........

i

r', +~ ',;~r *I~=o

I -",.l,,,i; I t !i I

11/t

",.,.",

Ill

POOLS Pc. Pi,Pn,P(]

ENZYMES Ep, E,[E P,]

Figure 1.

~, ~ "\',

(

J I...

2--d

I

I

I. i§

I

I

I

I "T::..12

!

" ...........

J

MESSENGERS Np,N,C, B,B~ Mp,M

..

.,,,"

i

,,~

I

GENES GCPI,Gp,G s

Flow chart showing d y n a m i c s of a cell nucleus (see Table 1)

GE represents a group of genes required for protein synthesis, Gp is used for the synthesis of the enzyme Ep (which is the polymerase for messenger M), gene Gs produces RNA fraction of the ribosome, while the protein fraction of ribosome is provided by the total protein E. Gene Gc produces transport RNA. The behavior of the cell system can briefly be described as follows. The complexing of messenger M with ribosome B yields template N; messenger M is produced by the interaction of the complex [EpPn] with gene GE. Enzyme synthesis results when templates interact with transport RNA and the amino acid pool complex [CPa]. The total protein E is used to convert internal pool P~ into RNA precursors and amino acids, and to convert external pool Pe into internal pool P~. Gene GE is the only gene where polymerase is operating in the messenger synthesis. The formulation of polymerase Ep follows the same line represented by total protein E synthesis. The external nutrient pool Pe is converted into an internal pool P~, where part of the pool P~ leaks out from the cell (via k3o), but the product is not associated with the external pool. The

SIMULATION OF CELL BEHAVIOR

431

principal regulating element in the system is the complexing total of protein with the internal pool (Eq. 22, Table I); this complexing is reversible and the ratio between the two rate constants k22, k32 determines the degree of regulation at this level. Additional regulating features are the complexing of polymerase with transport RNA (Eq. 21, Table I) and the complexing of polymerase with ribosomes (Eq. 23, Table I). 2.2 Mathematical Model of Cell. The following equations describe the mathematical model of the cell assumed in the paper; they are directly obtained from the chemical reaction equations of Table I using standard mass balance techniques: RNA Fraction of Ribosome dB' dt = k3GBPN -- k 2 B ' E

(1)

Ribosome dB d---[ = k 2 B ' E - k T B M - k s B M p - k14B + kloN[CPa]

-~- kllNp[CPa] - k 2 3 E r B + k33[EpB]

(2)

Transport RNA dC

d-"T = k4GcPn - k~CP~ + kloN[CP~] + kllNp[CP~] - k15C

- k21E~C + k31[EpC]

(3)

General Intercellular Metabolic Pool dP~ = l c l s E p e _ klgEP~ - k2oEP~ + k32[EPi] + k14B - k22EP~ dt

(4) + kl~C + k16 E + l~17E~ + k28 N + k29Np - k3oP ~

Polymerase-Transport RNA Complex d[EpC] = k21EpC - k3~[EpC] dt

(5)

Polymerase--RNA Complex dIE,B] dt

= lc2sE~B - kas[E~B ]

(6)

432

E.J. DAVISON

RNA Polymerase dEv = kllNv[CP.]

]~17Ev -- ]c21ErC + k31[EpC] --

-

dt

k23EvB

(7) + k33[EvB] -

blEvP n + koOE[ErPn]

Nucleotide Pool dP.

~t

= klgEP~ -

klEvP.

+ kmM

+ k13M v -

k4GcP~

(8) -

ksGpP~ -

k3GsP,

Polymerase-Nucleotide Complex d[EvP~! = klEvP ~ dt

(9)

koG~[EvP.]

Messenger RNA dM dt

= keGE[EvP.] -

k7BM

-

k12M + kloN[OPa]

(10)

Template dN - - ~ "= k 7 B M

-

kloN[CP~] -

(11)

]C2aN

Amino Acid Pool dPa dt = k 2 ~

-

(12)

kgCPa

Transport RNA-Amino Acid Complex d[CPa] =

kgCP~ -

Ir

-

(13)

kllNv[CP~]

dt

Total Protein dE

d--[ ffi k ~ ~

-

k~6E -

k22EPt + ka2[EPi] -

k2EB'

(14)

Protein-Intercellular Metabolic Complex d[EPi] dt RNA

= k22EP t -

(la)

ka2[EPJ

P o l y m e r a s e Cycle

Messenger RNA dMv dt

- -

= kaGvP n -

kaBM p + kalNv[OPa] -

k13M v

(16)

SIMULATION OF CELL BEHAVIOR

433

Template dNp = ksBMr dt State equations o f the cell.

- Ic11Np[CP~] - 1~29N~

(17)

For the sake of conciseness, the above equations

will be described by: i = f(x, u, k),

x(0) = x ~

0 < t _< t*,

(18)

a set of 17 simultaneous differential equations of the quadratic type, where x e R 17 is the state of the cell given by: xl = E ~ x2 = P~ x 3 = [EpP~]

x7 = M x8 = M s

x13 = E x~4 = P~

x9 = N

x l s = [EpC]

x4 = B ' x~ = B

Xao = N p x l l = P~

x16 = [EpB] x17 = [EP~]

x6 = C

x12 = [CPa]

(19)

where u e R 5 are the external inputs to the cell given by: ul = Pe u2 = Gs

ua = Gc u4 = Gs

u5 = G~ (20)

and which m u s t be determined, where k 9 R 29 is a vector of unknown rate constants for the cell which must be determined, and x(0) 9 R 17 are the initial conditions of the cell which m u s t be determined, t* is the division time of the cell, at which time cell division takes place, and must be determined. The above equations (18) describe the behavior of a single cell from the time of birth to the time at which cell division takes place (t* seconds later). 2.3 M a t h e m a t i c a l M o d e l o f a Growing Cell. To describe the behavior of a growing cell, it is necessary to introduce some t y p e of mitosis hypothesis (since the biological details as to when a cell divides axe not weU understood.). Consider the system (18); the following division hypothesis is made: D i v i s i o n hypothesis. L e t t* 9 [0, oo) be the time which globally minimizes the criterion function g(t) ~- IIx(t) - 2x(0)l I. Then if g(t*) < IIx(0)ll, cell division takes place at t = t*; ifg(t*) = IIx(0)l[, cell division does not take place. I f cell division takes place, the initial condition of the new generated cell is equal to 89 Figure 2 illustrates this division hypothesis. The motivation for such a hypothesis is t h a t a cell in steady-state behavior, i.e. with identical dynamic behavior existing from one generation of cell growth to another, will have the property t h a t J ( t * ) -- 0 < I]x(0)II, and x(0) = 89 L e t xt, t*, i = 1, 2, 3 , . . . denote the state and division time respectively of

434

E. J. DAVISON

Itx(t)- Zx(o)ll II.x(o)lt 0

t~

t

(ol

Jlx.(t) - 2x(o)lJ

II_x(o)1 ~

t

0

(b}

F i g u r e 2. Criterion for~deeiding w h e n m i t o s i s t a k e s place. (a) Cell divides a t t i m e t* in t h i s case. (b) Cell does n o t divide in t h i s case

the parent cell, 1st daughter cell, 2nd daughter cell etc. Then the behavior of a growing cell is t h u s described b y the following sequence of differentiM equations: f l = f(xl, u , k ) , i~+ 1 = f(xt+ 1 , u , k ) ,

xl(O) = x ~

x~+l(0 ) = 89

0 < t < t*, 0 < t < t*+l,

(21) i = 1,2,3 ....

where t*, i = 1, 2, 3 . . . . is found so t h a t min t,e[0, ~)llx,(t,) - 2x,(0)ll occurs at t~ = t*, assuming t h a t the constraint:

Ilx,(t *) - 2x,(o)ll < IIx,(0)ll, i = 1, 2, 3 . . . . (22) is satisfied. I f (22) is n o t satisfied for some value of i, say I, cell division will not t a k e place for i > I, and the cell stops dividing. Steady-state behavior in the cell is said to have taken place if there exists a scalar t* > 0 such t h a t : lim t* = t* t~ o~

lim IIx,+l(t) - x~(t)rl -- 0,

(23)

vte[O, ts*]

i--} ao

t* is called the steady-state division time of the cell in this case. I t will be assumed in the paper t h a t condition (23) m u s t be satisfied for a normal growing cell or for a cell which survives the effect of an input disturbance applied to it. I t will also be assumed, for the sake of simplicity, t h a t a cell which stops dividing, cannot survive, i.e. if (22) is not satisfied for all i, then the cell dies.* 2.4. Simulation of Cell Growth. I f it is assumed t h a t the parameters k, u, x ~ are known, t h e n the mathematical model of a growing cell (21) is completely * Note t h a t the topic of cell death is a very complex subject (Mitchison, 1965) and the fact t h a t a cell stops dividing does not, of course, necessarily imply t h a t it is dead, e.g. the cell may differentiate.

SIMULATION OF CELL BEHAVIOR

435

specified, and it is only a matter of numerically integrating the equations (18) and carrying out the 1-dimensional minimization search to obtain simulations of cell growth. In order to integrate the equations (21) an analogue or hybrid computer was not used becauseof scaling problems, and because there exist very successful algorithms for rapidly integrating large sets of differential equations (so called 'stiff' methods of numerical integration) on a digital computer, e,g. see Davison (1973). Since the differential equations (18) describing the cell were not 'stiff' in general, an adaptive 4th order R u n g e - K u t t a algorithm (IBM SSP R K , called RKGS) was used to integrate the equations (21). The following parameters were used in the algorithm: max step size = 0.05 see, max error = 10- 3 per unit step. The solution to (21) was calculated every 0.05 see for i = 1, 2, . . . , 15 with t ~ [0, T], T = 100 or until maxt~Eo, rl [Ix(t)[[ > 101~ the last precaution being taken to prevent overflow problems in the digital computer occurring. (Note that the growth behavior of (18), a set of unstable quadratic differential equations, is faster than an exponential, and so the numerical solution of (18): suddenly becomes extremely large.) Typically, the equations were integrated for t e [0, 20] using this last norm as a criterion. Using the above error criterion implied at least 3 significant figure accuracy in the computed solution which is more than adequate. This was verified by using a smaller error tolerance e.g. max error = 10 -6 per unit step. (A more accurate solution could have been obtained, but this would have been at the expense of more computing time.) Typical computation times required to integrate (18) on an IBM 370/165 digital computer for 0 < t < 20 were 0.68 sec real time. The 1-dimensional global minimization was carried out by an updating procedure.t 2.5. Identification of Parameters of Cell Model. At this stage in order to find a normal cell model, it is necessary to identify parameters k, u, x ~ (which are not necessarily unique) so that conditions (21) to (23) are satisfied, subject to the following constraints: k > 0,

u > 0,

xi(t ) > 0,

Vte[0, t*],

i = 1,2,3,...

(24)

where ~ > 0 means t h a t all elements of the vector ~ are constrained to be positive. (The last condition of (24) is necessary since all states of the cell are assumed to be positive.) This was done by using a Monte-Carlo search procedure for the parameters k > 0, u > 0, x ~ > 0 and then substituting these parameters into (18), and carrying out the cell growth process described by (21), (22) until condition (23) was satisfied. t U s i n g this u p d a t i n g proeodure, t h e division t i m e t* was f o u n d to a n a c c u r a c y of • 0.05 see; t h i s implies t h a t if t* is calculated to be e q u a l to only 0.05, t h e n condition (22) is n o t satisfied a n d t h e cell dies.

436

E.J.

DAVISON

In this identification process it was possible to take advantage of the inherent recovery property of the cell itself with respect to the initial condition pa~ameters x~ i.e. suppose k*, u* are the correct values to be used in the cell model and suppose t h a t x ~ is only approximately equal to the correct value given by x ~ say. Then, if the system (21)-(23) really does represent a real growing cell, one can think of x ~ ~ x ~ as resulting from a disturbance being applied to the cell. In this case the cell should 'recover' after a sufficient number of generations have elapsed, if the disturbance is not too severe (i.e. if x ~ is 'close enough' to x~ This means t h a t the initial conditions of the daughter cells will eventually approach the correct value x ~ after a sufficient number of generations have elapsed. This is illustrated in the following table (Table II) where TABLE II Table Showing How Initial Conditions of Cell States Were Identified Parent Cell

Initial Conditions Xl(0) Assumed (Listed from x I t o x17 ) 0.0170

1.36

0.384

0.218

0.304

0.114

0.334

0.244

0.078

0.066

1.14

1.03

0.240

0.850

0.0114

0.0197

0.332

1st Daughter C e l l

Initial Conditions x2(0) Obtained (Listed from x I to x17)

w

(t 1 = 5.75)

0.176

1.31

0.377

0.227

0.318

0,122

0.320

0.298

0.075

0.068

1.06

1.07

0.234

0,855

0.0125

0.0212

0.320 2nd D a u g h t e r C e l l (t 2 = 5.25)

Initial

C o n d i t i o n s x 3 ( 0 ) O b t a i n e d ( L i s t e d from x I t o XlT)

0.171

1.37

0.388

0.216

0.305

0.117

0.337

0.243

0.080

0.066

1.13

1.02

0.242

0,860

0.0118

0.0201

0.337 9th Daughter Cell (t 9 = 5.25)

Initial Conditions Xl0(0) Obtained (Listed from x I to x17) 0.170

1.37

0.385

0.218

0.305

0.115

0.335

0.245

0.078

0.066

1.15

1.03

0.241

0.850

0.0115

0.0198

0.333

SIMULATION OF CELL BEHAVIOR

437

k, u, x ~ have been chosen to be the correct values to be used in the model, except for the 1st element of x ~ which is chosen to be 1/10th of its correct value. In this oa~e, correct steady-state behavior (to 3 significant figure accuracy) is obtained after the 9th generation, with the initial conditions now being all correct. 2.6. Resul~ Obtained for Normal Cell. Using the identification technique presented in the previous section, the following results were obtained for the parameters of the cell model (see Tables I I I and IV): TABLE III Initial Conditions Obtained for Normal Cell (t~ = 5.25) Ep

[EpPn] B'

Pn

B

C

Mp

H

N

Np

Pa

[CPa]

E

Xl(O) x2(O) x3CO) x4(O) x5(O) x6(O) XTCO)Xs(O) XgCO)Xlo(O) Xll(O) x12(0) x13(O] 0.170 1.36 0.384 0.218 0.304 0.114! 0.334 0.244 0.078 0.066 1.14 1.03 0.240 [EpC] [EpB] [EPi]

Pi

x14(0) x15(0) x16(O) XlT(O) 0.850

0:0114 0.0197 0.332

TABLE IV Rate Constants and External Inputs Obtained for Normal Cell (t* = 5.25) k4

kS

k6

1,270

2. I00

0.847

10.93

k8

k9

klo

k12

1,053

3.408

3.833

kll 1.018

0.973

k16

k17

k18

I.030

0.109

1.606

k22

k23

k28 0.017

k1

k2

1.125

1.768

k7 3.178

k3

k13

k14

0.263

0.245

k19

k20

2.518

2.628

O.695

0.973

0.295

k30

k31 2.091

k32

k33

GH

0.880

1.278

0.i00

k29

0.054

0.500

GC

GB

0.100

0.100

kI5 1.077

k21

G

P O.I00

P

0

5.75

438

E.J.

DAVISON

Using these parameters the following results were obtained for the simulated cell model of (21), (22) (see Table V). TABLE V TableSho~ng Normal ~llReproduction Parent Cell

Initial

Conditions of Cell States

( L i s t e d from Xl(O) t o ~ 1 7 ( 0 ) ) .

0.170

1.36

0.384

0.218

0.504

0.114

0.534

0.244

0.078

0.066

1.14

1.03

0.240

0.850

0.0114

0.0197

0.332 F i n a l C o n d i t i o n s o f C e l l S t a t e s a t t 1 = 5.25 ( L i s t e d from x 1 (tl)

t o X l 7 ( t l ) ).

0.341

2.73

0.771

0.438

0.609

0.251

0.671

0.490

0.158

0.133

2.29

2.07

0.480

1.71

0.0251

0.0397

0.666 1st Daughter Cell

Final Conditions of Cell States at t 2 = 5.25 (Listed from

x 1 (t2) t o x 1 7 ( t 2 ) ) . 0.341

2.73

0.771

0.437

0.610

0.231

0.670

0.489

0.158

0.135

2.29

2.06

0.480

1.71

0.0231

0.0597

0.666

It can be seen that condition (23) is satisfied (to 3 significant figure accuracy) and that all states at time t* = 5.25 sec are equal to twice their initial value (to an accuracy of 1%). A plot of the transient behavior of the normal cell model for all states is given in Figures 3-9. It should be noted that the responses obtained agree qualitatively with Heinmet's results (Heinmet, 1966), but differ substantially in detail. 2.7. Examples Showing Normal Cell Dying. The following examples show the effect of changing some parameters of the cell model (from their nominal values) on the behavior of the cell system. Such changes of the model parameters

SIMULATION OF CELL B E H A V I O R

2.0

1.5

t I.o

0

F i g u r e 3.

!

3

5

t(sec)

?

V a r i a t i o n of P~, Pa, P~ for n o r m a l cell

2.0

1.5

[.0

t(ser

F i g u r e 4.

V a r i a t i o n of M , M s for n o r m a l cell

439

440

E. J. DAVISON

2.0

1,5

o

E

1.0

o

t(ser

Variation of E, Ep for normal cell

Figure 5.

2.0 -

B'

z~

t.5-

1.0

I

Figure 6.

3

5

t(see)

7

Variation of C, B, B' for normal cell

SIMULATION OF CELL B E H A V I O R

2.0

m

1.5

N E o z

1.0

tiser

F i g u r e 7. V a r i a t i o n of Np, N for n o r m a l cell

2.0

1.5 w N a E o Z

1.0

o

F i g u r e 8.

t($,r

V a r i a t i o n of [CPa],

[EpC] for

n o r m a l cell

441

442

E . J . DAVISON

2.0

EPi} ,~ 1.5 (Ep!

-g ,N E z 1.0

,B)

0

t(secl

Figure 9. V a r i a t i o n of [EPI], [E~P.], [EpB] for n o r m a l cell

correspond to applying various disturbances to the cell (see Tables VI and VII). In this case (Table VI) when the gene Ge is reduced by a factor of 10, cell division ceases at the 4th generation of cells and the cell eventually dies. Effect of

G e ---> ~ G ~

Division t Parent t .I

Cell

TABLE VI on Initial Conditions of Succeeding Generations of Cells (Cell Dies)

Time

Initial Conditions of Cell Xl(0)

x2(0)

x3(0)

x4(0)

...

x16(0)

x17(0)

0.170

1.36

0.384

0.218

...

0.0197

0.332

= 6.60

0.146

1.92

0.482

0.312

...

0.0174

I 0.332

r * = 7.40

0.145

1.64

0.440

0.389

...

0.0137

0 .237

t 3. = 7.20

0.132

0.840

0.223

0.4S6

...

0.0073s

0.075s

t 4 = 0.05

0.069

0.418

0.109

0.229

...

0 . 0 0 3 6 0 10.0370

* = O.OS t5

0.037

0.209

0.054

0.116

....

0.000171 0.0180

-

..,

Cell Dies

-

-

S I M U L A T I O N O F CELL B E H A V I O R

443

TABLE VII Effect of Pe -> ~ P e on Initial Conditions of Succeeding Generations of Cells (Cell Dies) Division

Time

Initial Conditions of Cell

Xl(O)

x2(O)

x3(O)

x4(O)

...

x16(0)

x17(0)

Parent Cell

0.170

1.36

0.384

0.218

...

0.0197

0.332

t I = 0.05

0.090

0.680

0.185

0.Ii0

...

0.00985

0.163

t 2 = 0.05

0.0490

0.340

0.090

0.055

...

0.00460

0.0795

t

Cell Dies

In this case (Table VII) the effect of reducing the concentration of the extracellular nutrient pool by a factor of 10, causes the cell to lose ability to divide even at the first generation. This corresponds to the well known phenomena in cellular biology, where it is known that cells are capable of growing only ahoy6 a certain minimal concentration of nutrient medium. 3. EFFECT OF Mg 2+ STARVATIONON THE CELL In order to verify the validity of the cell model, the following experiment was carried out: Let the ribosome content of the cell initially be put equal to zero (i.e. B = 0 at t = 0, or alternatively xs(0 ) = 0). Then how does the ribosome content of the cell vary with respect to rate of protein synthesis during recovery of the cell ? The experiment corresponds to an experiment McCarthy (1962) carried out in which cells of E. coli were starved of Mg 2 + content so that their ribosome content almost entirely disappeared. On restoration of Mg 2+, the ribosome content slowly returned to its normal value and the cells survived. The following results were obtained (see Table VIII) on simulating the system (21), (22) with all cell parameters equal to their normal values, except for x~(0) which is xs(0 ) = 0. It is seen that the cell recovers from ribosome deprivation (reaching steadystate behavior in 3 generations), thereby confirming McCarthy's results. On simulating the system (18) with all cell parameters equal to their normal values except for xs(0 ) which is xs(0 ) = 0, and on plotting B/0.313 x 100% vs E for 0 < t < t*, t* = 8.1, the following results were o b t ~ n e d (see Figure 10). In this case, it can be seen that the shape of the simulated cell recovery curve is almost identical to the shape of the experimental recovery curve of

444

E. 5. D A V I S O N

TABLE VIII Table Showing Steady-state Recovery of Cell Deprived of Ribosome Cell Division Time

Initial

Conditions of Cdlls

Xl(O)

xo(O)

x3(O)

x4(O) x5(O)

9. .

x16(0)

x17(0)

Parent Cell

0.170

1.36

0,384

0.218

...

0.0197

0.332

t 1 = 8.10

0.168

1.27

0.398

0.215 0.309

...

0.0230

0. 321

t 2 = 5.25

0.168

1.44

0.402

0.215 0.313

...

0.0202

0.356

t~ = S.O0

0.169

1.43

0.400

0.218 0.313

...

0.0200

0.352

t 4 = 5.00

0.169

1.43

0.400

0.217 0.313

...

0.0201

0.351

t

0

McCarthy's (1962). The only difference between the two curves is a bias effect, which can be accounted for b y the fact that the mathematical model of the cell is not intended to be a model of a specific cell, e.g.E, colt but only of a representative cell.

80

% Ribosome Content

6o

s lj

40

s, I

/ 20 /

O

i

i

[

]

O,2

0.4

0.6

0.8

_

Rote of Protein Synthesis (arbitrary traits)

F i g u r e 10. C o m p a r i s o n o f e x p e r i m e n t a l a n d t h e o r e t i c a l r e c o v e r y o f cell f r o m Mg 2 + s t a r v a t i o n . (a) C o m p u t e d c u r v e o b t a i n e d f r o m cell m o d e l . (b) E x p e r i m e n t a l c u r v e 4. E F F E C T

OF IONIZING RADIATION

ON A CELL

It is well known that if a cell is subject to a sufficiently large dose of ionizing radiation, the cell may eventually die. It is the purpose of this section to compare the simulated cell's response to ionizing radiation with experimental

SIMULATION OF CELL B E H A V I O R

445

results for the case of single-dose and split-dose ionizing radiation disturbances. In order to do so, it is necessary to relate the effect of ionizing radiation on a cell to the states of the cell itself. Since the phenomenon is biologically not completely understood, this will be done b y introducing a radiation effect hypothesis.

Radiation effect hypothesis. A dose D of ionizing radiation at time t affects a cell b y destroying a certain fraction of all the components (states) of the cell at time t; the larger the dose D is, the larger is the fraction of the cell's components (states) destroyed at time t. The motivation for the hypothesis is that experimental evidence has indicated that ionizing radiation actually destroys part of the various components of the cell being radiated (Little, 1968). It will be assumed then that the effect of ionizing radiation on a cell at time t = 0 is to cause the initial conditions of all states of the cell to be multiplied by a number r, 0 < r _< 1 called the radiation factor; the more severe the dose of radiation, the smaller is the radiation factor. 4.1 Single-Dose Experiments. Assume now that a single dose of ionizing radiation of various degrees of intensity is directed on a cell at time t = 0. Then the results of simulating the system (21), (22) using the normal parameters of the cell, except for the initial states which are all multiplied by a number r, the radiation factor, are given in Tables I X and X. TABLE IX Ionization-Radiation Ex ~eriment (Single Dose Given) Radiation f a c t o r r

D i v i s i o n Time of R e s u l t i n g S t e a d y - s t a t e Cell Growth

0.i

cell dies

0.2

cell dies

0.60

cell dies

0.65

9.65

0.70

8.60

0.75

8.10

O. 80

7.10

0.85

6.65

0.90

6.30

0.95

5.75

1.00 normal initial conditions

S.25 normal c e l l

446

E. J. DAVISON TABLE X Some Typical Resu]ts Obtained for Single Dose of Ionizing Radiation Experiment Radiation Factor

Division

Time

Initial

C o n d i t i o n s o f C e l l S t a t e ( L i s t e d in

o r d e r Xl(O ) to x17{0))

r - 0.2

r - 0.6

t 1 : 0.05

!Cell

t~ o 0.05

t

t [ = 14.1

0.176 0.520 0.166 0.240 0.185 0.072 0.164

d o e s not d i v i d e .

Cell dies.

0.032 0.027 0.510 0.750 0.084 0.70

0.116

0.0081 0.0144

0.107

t 2 : O.OS Cell does not divide.

Cell dies.

t~ : 0.05 r : 0.65

t~ : 14.1

0.182 0.790 0 . 2 5 3 0.225 0.239 0.094 0.234 0.159 0.052 0.041 0.730 0.850 0.140 0.790 0.011 0.019 0.187

t~ : I0.1 t~ ffi9.95 t; = 9.85

t~ : 9.75

Steady-state Behaviour Reached

{

t 6 : 9.65

t7

9.65

0.177 0.840 0.260 0.228 0.239 0.092 0.240 0.173 0.052 0.044 0.790 0.870 0.140 0,790 0.010 0,017

0,193 t~ = 0,10

r =0.75

0.180 0.920 0.285 0.229 0.253 0.098 0,250

0,189

0.058 0.049 0.835 0.905 0.160 0.805 0.0108 0.0184 ).215 Steady-state Behaviour Reached

t

t 2 = 8.10 t5

8.10

0.175

0.945

0.286 0.226 0.251 0.096 0.261

0.189

0.0585 0.0489 0.865 0.895 0.162 0.805 0.0102 0.0171 O. 220

I t is seen t h a t if the radiation dose is too severe (i.e. for r ___ 0.60) t h e cell dies a n d for w e a k e r radiation doses (i.e. for r > 0.65} t h e cell survives w h i c h agrees w i t h e x p e r i m e n t a l results (Little, 1968). It m a y be observed in Table X t h a t for r = 0.6, the cell divides once before it dies, whereas for r = 0.2, t h e cell dies

SIMULATION OF CELL BEHAVIOR

447

immediately without dividing. This observation agrees with the following remark taken from Little (1968): "A rather perplexing feature that remains however, is the ability of lethally radiated cells to undergo several normal division cycles before mitosis is inhibited." 4.2. Split-Dose Experiments. Consider now the effect of applying only 89of the radiation dose at t = 0 and the other 89of the dose A sec later. Experimentally it has been found t h a t the cell in this case is much less affected by the radiation as compared to a single dose of radiation. See for example Elkind et al. (1964), Sinclair and Morton (1964), who study the effect of X-ray radiation on Chinese Hamster cells. r I = 0.75 r2= 0 . 7 5 division time 9 5 . 2 eec

.

.

of split oos,

9 F

9.! -.

:

.

6.6

~

6,6

9

/ 2 5 r(,.,: ;.".///,)..'//A, ,I~2.;(~,'.;:;.....;~,,J//~'? o

~

~

s

A T i m e Intervo|

s Between

,o

,2 Doses

,4

,s

,s

2o

(See)

Figure 1l. Computed split dose recovery function for a single cell subject to ionizing radiation obtained from cell model. Survival is plotted as a function of the time interval between the two doses rz, r2. Normal steady-state division time of the cell is 5.2 see I t is of interest to see if the cell model considered also behaves in this way. The following results were obtained on simulating the cell model (21), (22) assuming that all cell parameters are normal, except for the initial conditions which are all multiplied by the radiation factor rz = 0.75 and, then A sec later, multiplying all resultant states of the cell by the radiation factor r2 = 0.75 (see Figure 11). This corresponds to applying a total dose with a radiation factor of (0.75) 2 = 0.563. This can be compared with the experimental results obtained by Sinclair and Morton (1964) given in Figure 12. As can be seen, the computed response curve is very similar to the experimental curve; in particular, it has a minimum occurring at a time somewhat larger than the division time, which agrees very well with the experimental results. (Note t h a t the recovery curve of Figure 11 is only for one cell; for a large number of cells, taking into account the statistical variation of the cells, the qualitative predicted curve of Figure 12 will occur.) 4.3. Mechanism of the Split-Dose Recovery Curve. The following explanation is given of the shape of the curve given in Figure 11. For A small (0 < A < 2.5),

448

E.

J.

DAVISON

the two doses of radiation are very close together, and so the overall effect of radiation approximately corresponds to a radiation factor of 0.75 • 0.75 -0.56, which from the single dose experiment (Table IX) is known to cause the cell to die. Therefore for A small enough (i.e. 0 < A < 2.5) the cell dies. As A gets larger (2.5 < A < 9.1), the effect of the 2nd dose on the cell takes place at a much later time ~ > 2.5, and at this time ~ the natural recovery (growing) properties of the cell have started to minimize the effects of the 1st dose on the cell; in particular, the magnitude of the states of the cell are much 0.05 DI = 730 rods D2=750 rods division time = 5.5 hrs.

0.04 0.03 Fraction Surviving Split Dose 0 0 2 Radiation "

/ ,~/,,.~.

0.01

/

~"~..~ v I

O

I

~

%

~

/~r

[

-- ~ I

. . , "Experimental I

1

,," r ,r \ "~

~

I

2 4 6 8 1O t2 14 16 A Time Interval Between Doses (hrs)

predicted Resutt (Qualitative Only) I

18

9

20

Figure 12. Experimental split dose recovery function in cells subject to ionizing radiation (Sinclair and Morton, 1964). Survival is plotted as a function of the time interval between the two doses D1, D2. Normal steady.state division time of the cell is 5.5 hr larger than they were for t < ~ and are approximately equal to the initial conditions of a normal cell for some t e [2.5, 9.1]. (See Table X: recall t h a t the terminal state of the parent cell is equal to twice the initial conditions of the daughter cell given in Table X.) Therefore when the 2nd dose of radiation occurs at time ~, it is approximately equivalent to just a single dose of radiation on a normal cell with a radiation factor of 0.75, which, from the single dose experiment (Table IX), is known to be mild enough to not affect the survival of the cell. Therefore, for A large enough (i.e. 2.5 < A < 9.1), the cell will survive. As A gets larger yet (9.1 < A < 10.6), the cell will eventually divide with only the 1st dose of radiation having affected it, and thus a new daughter cell appears at time 9.1 (see Table X). I f a 2nd dose of radiation is now applied very shortly after this division time (in particular if 9.1 < A < 10.6), then the new daughter cell will receive a dose of radiation with a radiation factor of 0.75; however, at this time (from 0 to 1.5 sec after the division time at 9.1), the new daughter cell is in a much more fragile condition than the corresponding normal cell would be (in fact, the magnitude of the states of the new cell at this time are only approximately ~ of their normal values--see Table X for the initial con-

SIMULATION OF CELL BEHAVIOR

449

ditions of the daughter cell), and as a result, it cannot survive the 2nd dose of radiation. Thus for A slightly larger than the time at which cell division takes place (i.e. 9.1 < A < 10.6), the cell dies. As A gets even larger (10.6 < A < 17.2), the natural growth properties of the daughter cell have started to minimize the effects of the 1st dose on the cell, and eventually the daughter cell becomes 'strong enough' (i.e. the states of the daughter cell becoming almost normal) to survive a 2nd dose of radiation occurring after 10.6 sec. Thus for A large enough (10.6 < A < 17.2), the cell will survive. As A gets larger yet (17.2 < A < 18.2), the cell will eventually divide and thus a second daughter cell appears at time 17.2 (see Table X). If a 2nd dose of radiation is now applied very shortly after this time (in particular, if 17.2 < A < 18.2), the new daughter cell will receive a dose of radiation and since the new daughter cell is in a more fragile condition1: than the normal cell, it cannot survive this 2nd dose of radiation. Thus for A slightly larger than the time at which cell division takes place the second time (17.2 < A < 18.2) the cell dies. For A extremely large (h > 18.2), the natural growth properties of the 2nd daughter cell have now minimized the effects of the first dose of radiation and now the daughter cell is 'strong enough' {i.e. the states of the daughter cell becoming almost normal) to survive a 2nd dose of radiation occurring after 18.2 sec. Thus for A extremely large (A > 18.2) the cell will survive. Remark. It can be seen that an explanation as to the shape of the split-dose recovery curve can be made without the necessity of introducing extra hypotheses at this stage, i.e. the complex recovery curve is due simply to the way the biochemical entities of the cell interact; one does not have to introduce various 'healing' hypotheses (Little, 1968) to explain the curve. The analysis indicates that the first local minimum appearing in the recovery curve will occur at a time somewhat larger than the division time of the cell. 5. EFFECT OF DISTURBANCES ON THE CELL: ' T u M o R - C E L L ' BEHAVIOR OF THE CELL

This section deals with the systematic study of applying various disturbances to the cell and then observing the effect of these disturbances on the cell's behavior. In this case, disturbances can affect the cell by changing any of the rate constants, the external inputs, or the initial conditions of the cell. The following results were obtained on simulating the system (21), (22) on using the normal $ Note t h a t this is not always true; it depends on the relative doses rl, r2 whether the 2nd minimum in the survival curve appears.

450

E.J. DAVISON

parameters of the cell, except for one parameter only, which was changed by a factor of 10 or ~ (see Table XI). In this case it can be seen that for some disturbances, the cell dies, for other TABLE XI S u m m a r y of Results O b t a i n e d when I n i t i a l Conditions o r P a r a m e t e r s o f Cell are Changed x I0

Parameter P e

Steady-state *

**

1.9

Steady-state Division time t s ! dies dies dies

dies *

dies

(1.9)

5.9 (3.6)

1.8 (1.8)

*

1.O ( i . 0 )

dies

*

0.75 (0.75)

4.1 (5.1)

dies

dies

dies

dies

k4 kS k6 k7

(3.9)

dies

Gb Gc Ge G P kI k2 k3

1.7

xO.l

Division Time =w'-'t s

6.0 (3.7)

1.8 ( 1 . 8 )

*

1.9 (1.9)

dies

*

1.1 (1.3)

dies

3.6 (2.1)

2.3 (1.8)

k8

k9

*

1.9

(2.1)

2.3 ( i . i )

k10 kll

*

1.6

(1.7)

7.0 (2.1)

*

1.6 (1.6)

5.8 (2.3)

k12 k13

dies

2.8 (3.2)

2.3 (2.4)

6.2 (4.6)

k14 kls k16 k17 k18 k19

dies

dies

2.2 6i.7)

S.l (3.9)

dies

2,8 (4.0)

*

k20 k21 k22 k23 k28 k29 k30 k31 k32 k35 Xl(O)

x2(o)

s.2 (5.3)

S.6 (4,3)

*

1.7

(3.8)

**

dies

dies dies

dies

5.3 (1.3)

1.2 (1.2)

12.4 (2.3)

1,7 (1.9)

S.S (1.3) S.6 (5.3) s.3 (s.4)

1.2 (1.2)

4.4 (5.0)

4.4 (s.1)

dies

S.2 (5.4)

S.6 (2:2)

S.S (2.2)

7.4 (2.7)

7.4 (2.7)

.5.8 (2.6)

5.8 (2.6)

dies 1.2

5.3 S.2

SIMULATION Table

x3(0): x4(0) xs(0) x6(0) xT[0) x8(0) xg(O) Xl0(0)

*

*

1.2 1.2

*

1.2

*

1.6

*

1.5 2.7 2.5

*

1.9

d ies

xls(0)

*

1.2

x16(0)

*

1.2

x17(0)

*

1.9

45I

5.0 5.2 5.0 5.3 5.i 5.3 5.2 5.3 5.3 5.2 $.3 5.2 5.3 5.2 5.2

1.9

*

BEHAVIOR

(Continued)

dies 2.7

Xn(0)

x12(0) x13(o) x14(0)

11

OF CELL

( ) Refers to the d i v i s i o n time achieved when the s t e a d y - s t a t e d i s t u r b e d c e l l has i t s perturbed *

parameter

changed back

to i t s normal v a l u e .

Refers to the fact c h a t c e l l i s behaving l i k e a f a s t growing t u m o r - c e l l ,

disturbances the cell basically remains normal, and for other disturbances yet, the cell's steady-state division rate increases (by as much as a factor of 7!). I t is also to be observed t h a t those disturbances (marked with * in Table XI) which cause the cell's steady-state division rate to significantly increase in Table XI are irreversible, i.e. once the disturbed cell has become tumor-like, there is no change in the behavior of the cell if the perturbed parameter of the cell is now changed back to its original normal value. (This statement is not quite true for the two disturbances marked with an ** in Table XI). A sensitivity study of some typical disturbances which cause the cell's stea~tystate division rate to significantly increase (disturbances marked with * in Table XI) is given in Table XII. As can be seen in Table XII, some disturbances (in particular, those marked * in Table XII) cause the cell to very suddenly increase its steady-state division rate. This is illustrated in more detail for the disturbance affecting xa(0) in Table XIII. Table XIV gives a comparison of the physical size of some of the typical fast growing disturbed cells of Table XI with a normal cell. As can be seen, the fast growing disturbed cells are much larger in magnitude than the normal cells (some states being 50 times larger than the corresponding normal state!). This observation is true for all the fast growing cells marked * in Table XI. I t can be seen that these fast growing cells are behaving in a way which is very

452

E. J . D A V I S O N

TABLE X I I Summary of Results for Sensitivity Study of Some Typical Fast Growing Tumor Cells Obtained when Initial Conditions or Parameters of Cell are Changed

Disturbance Factor

Steady-state Division Times (t*) for Various Initial Conditions or Parameters Pe

Ge

kl

k2

1.0

0.75

i.I

1,2

1.2

1.2

1.2

1.2

1.2

0.80

1.9

1.8

2.4

2.1

4~7

2.5

2.4 2.9

k7 x2(O) x6(O) x7(O) x8(O) Xl5(O) x16(O)

x 10

1.7

1.9

x5

2.1

2.7 0.80

x4

2.5

3.0

0.85

717

2.2

2.2

2.9

2.6

4.4

3.0

x3

2.8

3.6

3~0

6.5

2.6

2.9

3.6

3.4

4.3

3.5

3.6

x2

3.5

4.4

3.9

5.7

3.2

4.1

4.5

4.5

4.8

4.3

4.5

4.9

4.9

x 1.5

4.0

4.8

4.4

5.2

4.0

4.8

5.0

5.0

5.0

x I.I

4.9

5.1

5.1

5.3

4.9

5.2

5.1

5.2

5.2

5.1

5.2

xl

5.3

5.3

5.3

5.3

5.3

5.3

5,3

5.5

5.3

5.3

5.3

Disturbance Factor

Division Times its) G p

k21

*'Sudden change occurs at this point. xO.1

L8

1.2

xO.2

3.0

2.6

xO.4

4.4

3.9

xO.6

4.7

4.7

xO.8

5.0

5.1

xO.9

5.2

5.2

xO.9S

5.2

5.2

xl

5.3

5.3

TABLE XIII Sensitivity Study for Fast Growing Tumor Cell Resulting from Change of Initial Condition on xs(0) Disturbance Factor x I0 x 9 x 8 x 7 x 6.5 x 6.0 Steady-state Divi-

1.2 1.2 1.3 i : 5

1'.50

1.65

5.5 x 5 x 4 x 2 x 1.5 x 1 . L x 4.954.74.44.9

5.0

5.2

1

5.3

sion Time (t;) * Sudden change o c c u r s a t t h i s p o i n t .

similar to tumor cell behavior. Table XV shows the w a y these fast growing 'tumor' cells approach steady-state behavior in approximately 4 generations. Figure 13 shows how a typical fast growing 'tumor' cell behaves during the time interval 0 < t < 1.2, assuming it has the initial condition given in Table XIV.

SIMULATION OF CELL BEHAVIOR

O

~

~

o

~

o

o

o

o

o

~

o

9

o

~.~ r

o

~

o

o

o

o

~

~

o

~

o

~

o

o

o

.,-* o

n

"~

~ o

453

454

E. J. DAVISON 2 XI

x2

O

x3,

I

i

t

I

o

f

i

r

i

t

i

t

o

]

t

0

I

f

/ x4

x5

o

i

t

X7

X8

o

L

i

;/ o

Xll ;

I

t

XlO o

1

i

X12 0

I

x16

I

t

xt~ f

X14 I

t

I

,

l

t

Xl3

x9

0

1

0

I

t

XI'/' I O

1

t

F i g u r e !3. F i g u r e s s h o w i n g s t e a d y - s t a t e t r a n s i e n t b e h a v i o r o f f a s t g r o w i n g " t u m o r " cell r e s u l t i n g f r o m d i s t u r b a n c e xs(0)~-> xe(0) • 9. ( O r d i n a t e s are normalized)

5.1. Perturbations which cause Cell to Become Tumorous. The following is a summary of all parameter perturbations which cause the cell to suddenly change into a fast growing 'tumor-like' cell with the following properties: (i) the disturbed cell is physically larger than the normal cell, (ii) the division time of the disturbed cell is greatly speeded up, (iii) the energy requirements of the disturbed cell are much greater than those of the normal cell, (iv) the disturbed cell behaves in an irreversible way, i.e. it is not possible to apply an opposite parameter change to the disturbed cell so that it returns to its normal behavior. The cell becomes 'tumor-like' when any one of the following changes in the normal parameters of the cell occurs:

(1)

Pe --> Pe x 10

(2)

Ge ~

(3)

Gv ---->G p x ~

(4)

Pn(O)-+ P.(O) x 10

Ge x

10

SIMULATION OF CELL BEHAVIOR

(5)

[ErP.](0 ) --> [ErP.](0 ) x 10

(6)

c(o)-~c(o)

(7)

M(0) --> M(0) x 10

(8)

Mr(0 ) - ~ M r ( 0 ) • 10

(9)

N(0) --> N(0) x 10

(lO)

N d o ) - + N d o ) • lO

(11)

E(0) - ~ E(0) • 10

(12)

[E~C](O) ~ [ErC](0) • 10

(13)

[E,,B](0) --> [ErB](0 ) x 10

(14)

[EP,](O)-+[EP~](O) x 10

(15)

Ep + P~ ~'>[EpP~];

(16)

B' + E ~2~B;

(17)

G~ + P n k S > G , + Mr;

(18)

GE + [E~P.] k6> GE + E r + M;

(19)

B + M k~N;

(20)

C + Pa k">[CPa];

(21)

N +[CPa] k~O>B + C + M + E;

(22)

Np + [CPa] ~'~t> B + M r + C + Er;

(23)

Pe + E ka~>pi + E ;

(24)

E r +C.~2'

(25)

E + Pi k ~ [EPi];

(26)

E r + B ~ ' - [ErB];

k31

• lO

k~3

kl-->kl • 10

k2-->/c 2 x 10 k5-->/c5 • 1~ k6-->k s • 10

]c7--->k 7 • 10

[EpC];

R32

/~23

455

]r

• 10

klo--> klo • l0 kll ---->kll x l0

k18--+/c18 • l0 k 2 x ~ k 2 1 • 16 k22-->k22 • /c23--> k23 • ~ .

It would be interesting to determine those external disturbances in the cell which cause the above perturbations in the cell's parameters to occur--this then would give insight into what t y p e of external disturbances in a cell cause it to behave in a tumorous way. This is a future research area.

456

E.J.

DAVISON

6. CoNcLusions This paper has obtained a mathematical model of a dividing cell and studied various behavioral patterns of the resulting cell growth, when the cell is subject to different disturbances. The following main results have been obtained: (a) The cell model's behavior has agreed with experimental behavior for the case of a cell subject to Mg 2+ starvation and then allowed to recover to nominal values. (b) the cell model has agreed with experimental behavior for the case of a cell subject to a split-dose of ionizing radiation. In this case the cell model has explained a mechanism for the rather complicated survival versus time interval plot of cells subject to a split-dose of ionizing radiation. In explaining this mechanism, there is no need to introduce extra assumptions such as 'healing effects' occurring in the cell. (c) The cell model has predicted that fast growing 'tumor-like' behavior can occur in the cell when certain type of disturbances affect it; in such cases the physical size of the cell is much larger than its nominal value (typically b y a factor of 10), the division time of the cell is greatly speeded up (typically b y a factor of 5), and the energy requirements for the disturbed cell are much larger than for the nominal cell. Moreover the behavior of these 'tumor-like' cells is irreversible, i.e. once the disturbed cell has become tumor-like, there is no change in the behavior of the cell if the perturbed parameter of the cell is now changed back to its normal value. Since the mathematical model did not include any effects of extra-nuclear behavior in a cell (e.g. the cytoplasm) or interaction of cells with neighboring cells, this implies that tumor behavior (or malignant behavior) can be bought about b y disturbances occurring solely in the chemistry of the nucleus of the cell, and hence is a very fundamental property of a cell; i.e. one does not need to introduce any more complexities into the mathematical model of the cell in order to achieve tumor growth. It would obviously be of interest to verify experimentally these predictions. This is an area of future research. The most critical assumption made in the mathematical model of the cell is the 'division hypothesis' which states that the cell will divide at a time t* when the norm given b y [Ix(t*) - 2x(0)]] is minimized, provided ]Ix(t*) - 2x(0)]] < I]x(0)ll, where x(t)is the state of the cell. This hypothesis was made because the biological details as to when exactly a cell divides are little understood. The motivation for this hypothesis is that a cell in steady-state behavior (i.e. dividing at a constant rate with identical dynamic behavior occurring between different generations of cells) has the property that mint. ]Ix(t*) - 2x(0)l I = 0.

SIMULATION OF CELL BEHAVIOR

457

This hypothesis gives a model of a cell which agrees with experimental evidence for the case of a cell subject to split dose ionization radiation. Future work on the problem would involve changing this hypothesis and introducing new hypotheses which possibly have a more detailed biological basis for motivation. The above study has indicated that the mathematical study of cell processes is a method which can be used to obtain insight into the behavior of biological processes. It has the significant advantage that simulations of complex behavior can rapidly be carried out on a digital computer with all interior states of the process being monitored. For example, the total time required to simulate a cell dividing once was only 0.7 sec of real time on an IBM 370/165 digital computer; this includes the time to plot all 17 states of the cell as a function of time and perform data collection. This author wishes to thank E. Fabian, P. Wong, E. Mak and J. Phillips for assistance in the computer modelling of the cell system. This work has been supported by the National Research Council of Canada under Grant A4396. APPENDIX TERMINOLOGY AND SYMBOLS Pools Pe Pt Pa Pn

= = = =

Extracellular nutrient pool General intracellular metabolic pool Amino acid pool for protein synthesis Nucleotide pool for R N A synthesis

Enzymes E -- Total protein Ep = R N A polymerase for messenger R N A (M) synthesis Ge,r$e8

GE = Genes for messenger R N A (M) synthesis Gp = Genes for messenger R N A (Mp) synthesis Gs = Genes for the synthesis of R N A fraction of ribosome

Gc = Genes for transport R N A (C) synthesis

Messengers M = Messenger (RNA) for protein (E) synthesis Mp = Messenger (RNA) for Ep synthesis B" = R N A fraction of ribosome B = Ribosome C = Transport R N A N = Ribosome and messenger complex for protein (E) synthesis (template) Np = Ribosome and messenger complex for Ep synthesis (template)

Intermediate Complexes [CPa] [EpB] [EpPn] [EP~] [EpC] X

LITERATURE Davison, E. J. 1973. "An Algorithm for the Computer Simulation of Very Large Dynamic Systems." Automatica, 9, 665-675. Diamant, N. E., P. K. Rose and E. J. Davison. 1970. "Computer Simulation of Intestinal Slow-Wave Frequency Gradient." A m . J. Physiol., 219, 1681-1690.

458

E. J. DAVISON

Elkind, M. M., T. Alescic, R. W. Swain, W. B. Moses and H. Sutton. 1964. "Recovery of Hypoxic Mammalian Cells from Sub-Lethal X - R a y Damage." Nature, Lond., 202, 1190-I193. Glueck, A . R . 1969. "Simulation of Cell Behavior." Private communication. Heinmets, F. 1964a. "Analog Computer Analysis of a Model-System for the Induced Enzyme Synthesis." J. Theoret. Biol., 6, 60-75. 19641). "Elucidation of Induction and Repression Mechanisms in Enzyme Synthesis b y Analysis of Model System with the Analog Computer." I n Electronic Aspects of Biochemistry, pp. 413-479. New York: Academic Press. 1966. Analysis of Normal and Abnormal Cell Growth. New York: Plenum

Press. Himmelblau, D. M. and K. B. Bischoff. 1968. Process Analysis and Simultation. New York: Wiley. Jansson, B. 1968. "Mathematical Description of the Growth of Tumor Cell Populations with Different Ploidi-Compositions," 6th Annual Symposium in Math and Computer Science in the Life Sciences, Houston. Little, J . B . 1968. "Cellular Effects of Ionizing Radiation." New England J. ~ledicine, 278, 369-376. Mallette, M. F., C. P. Clagett, A. T. Phillips and R. L. McCarl. 1971. Introductory Biochemistry. Baltimore: Williams & Williams. McCarthy, B . J . 1962. "The Effect of Magnesium Starvation in the Ribosome Content of Escherichia Coll." Biochim Biophys Acts, 55, 880-888. Medoff, G. and M. N. Swartz. 1967. " D N A - - S t r u c t u r e and Enzymatic Synthesis." New England J. Medicine, 277, 728-737; 276, 788-795. Mesarovic, M. D. (ed.). 1968. Systems Theory and Biology. New York: Springer. Mitchison, J . M . 1965. The Biology of the Cell Cycle. Cambridge University Press. Palmby, F. V., E. J. Davison and J. Duffin. 1964. "The Simulation of Multi-Neurone Networks: Modelling of the Lateral Inhibition of the Eye and the Generation of Respiratory Rhythm. Bull. Math. Biology, 36, 77-89. Pollard, E. 1960. "Theoretical Aspects of the Effect of Ionizing Radiation on the Baterial Cell." Am. Naturalist, XCIV, 71-84. Sinclair, W. K. and B. A. Morton. 1964. "Survival and Recovery in X-irradiation Synchronised Chinese Hamster Cells." Cellular Radiation Biology Proc., 18th Ann. Symp. .Fund. Ca~vcer Res., Univ. Houston, Texas. Teir, H. and T. R y t 6 m a a (eds.). 1967. Control of Cellular Growth in Adult Organisms. New York: Academic Press. Warner, J. R. and R. Soeirs. 1967. "The Involvement of RNA in Protein Synthesis," New England J. Medicine, 276, 563-570; 276, 613-619; 276, 675-680. Weinberg, R. and ]3. P. Zeigler. 1969. "Computer Simulation of a Living Cell: Multilevel Control Systems." Report 08228-17-T. Ann Arbor, Michigan: University of Michigan. Yeisley, W. G. and E. C. Pollard. 1964. "An Analog Computer Study of Differential Equations Concerned with Baterial Cell Synthesis." J. Theoret. Biol., 7, 485-503. RECEIVED 1-17-75