- Email: [email protected]
Simulation of cell behavior: Normal and abnormal growth
B U L L E T I N OF MATHEMATICAL BIOLOGY
VOLUME
37, 1975
S I M U L A T I O N OF CELL B E H A V I O R : NOI~MAL A N D A B N O R M A L G R O W T H
•
E. J. DAVmON 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 disturbanees, in particular when (i) a Mg~+ 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 are perturbed. The cell model results obtained agree well with experimental results for the Mga + 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.
h INTRODUCTION The purpose of this p a p e r is to find the possible m e c h a n i s m which explains w h y a n o r m a l cell s u d d e n l y becomes 'tumor-like'. This will be done b y obtaining a m a t h e m a t i c a l m o d e l of a growing cell living in a culture medium, a n d t h e n s y s t e m a t i c a l l y e x a m i n i n g the b e h a v i o r of the cell w h e n various disturbances are applied to different parts of the m a t h e m a t i c a l model of the cell. The m o t i v a t i o n for a p p l y i n g such a model a p p r o a c h is t h a t m a t h e m a t i c a l s y s t e m techniques h a v e been highly successful in explaining a n d predicting the b e h a v i o r of technical-industrial processes (see e.g. H i m m e l b l a u and Bischoff, 427
428
E . ft. D A V I S O N
1968), and recently have given insight into the mechanism of various biological processes (see e.g. Mesarovic, 1968; Diamant et al., 1970; P a l m b y 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. I-Ieinmets (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 8e, 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 t r u e 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 case 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 ceils 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 ~+ 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. ~OI~MAL 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.
E p + Pn " ~ [EpPn]
2. B' + E k-~2B k~ 3. GB + Pn _~o GB + B' k4 4. GC + Pn .-:-4. Cp + ,'.Ip kS
5. Gp + Pn -~" Gp'+ H
P
6.
k6 ~1 GE + [EpPn] ~ G]] + Ep +
7.
k7 B + !,I ~ N P
i
!o.
B kI"~4P, 1
15. C kl-~S P. 1 16. E kl-~6 p. 1.
17. g
p
k 17p.
1
18. P + E kl-~8P. + E 19. P.l + E - ~ Pn + E 20. P. + g k2-~O P
5.. P
, 1
c+P#[cPa]
lO,'N + [CPa] ~
14.
B+ C+H+E
ii. Np + [CPa] k'~11B + ~,IP + C + EP
21, Ep
+
C
k31 [EpC]
k22 22. E + Pi ~32 [EPi] 23. Ep + Bkk~3 [EpB]
12. ~Ik-~ 12 P n 13. M k-~13 P p n
+ E a
24. N-~-~ P.
1
25. N - ~ P . p I 26. p. k'-~ X 1
430
E. J . D A V I S O N
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.
RNA Synlhes~s
k9
~i
/ "-- "i
1 I l
Pool PI leal~age
7i kl
///
~<
/
/
s
INTERMEDIATE COMPLEXES
r
/
"
/ Messenll~,
ti ~~,-+-ffl-t----~
Pi
o
.7!1 -I+L
1 '
\
I
/,,.+ ..-- ~
x I +'I'~ ~-
.....
.
.
; 1<' -.*'.I,=, ~Pj
.
.
.
.
~,
t£,,.
It
.
I
=
. . . .
I
/"
1
~',.-
I
. ¢ .......................
,% POOLS Pe, PI,PI~,~o
ENZYMES Ep,E,[EPI]
MESSENGERS
Np,N,C,BiB%Mp,M
GENES GOGI~,Gp,GE
Figure 1. Flow chart showing dynamics of a cell nucleus (see Table 1) G E represents a group of genes required for protein synthesis, Gp is used for the synthesis of the enzyme E , (which is the polymerase for messenger M), gene G B produces I~NA 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 G E. Enzyme synthesis results when templates interactwith transport I~bTA and the amino acid pool complex [CPa]. The total protein E is used to convert internal pool Pi into I~NA precursors and amino acids, and to convert external pool Pe into internal pool Pi. Gene G E 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 k30), but the product is not associated with the external pool. The
S I M U L A T I O N OF C E L L B E H A V I O R
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 k~, k3~ 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 M o d e l o f 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 Fr~ction of l~ibosome dB' d--T = lc3GBP~ - lc2B'E
(1)
Ribosome dB d---{ = k 2 B ' E -
k T B M - k s B M p - k l ¢ B + kloN[CPa]
+ k l l N p [ C P ~ ] - k23E~B + k33[EpB]
(2)
Transport RNA dC
d-T = 1QGcP~ -
kgCP~ + k l ° N [ C P a ] + kllNP[CPa] - k15C
_ lc21E~C + lc3~[E,C]
(3)
General Intercellular Metabolic 1)o01 d P i = k l s E P e - k 1 9 E P i - k2oEP~ + lc3e[EP~] + k14B - ]c22.EP ~ dt
(4) + k15C + k16E + k~TE~ + ]c28N + lc29Np - k3oP~
Polymerase-Transport RNA Complex dIE, C] = lc21E, C _ k81[EpC] dt
(5)
Polymerase-RNA Complex dIE,B] dt
lc23EpB - lc33[EpB ]
(6)
432
E.J.
DAVISON
gNA Polymerase dEn = knNp[CP~] dt
-
]QTEp -
+ k33[EpB] -
k 2 1 E ~ C + ]c3a[EpC] -
k23EpB
(7)
k l E p P ~ + k6GE[E~,Pn]
Nucleotide Pool dP,~ d t = ]clgEPi -
klEpP~
ksGpP n -
+ ]c12M + k 1 3 M p -
k4GcPn
(8)
k3GBP ~
Polymerase-Nucleotide Complex d[E,P~] dt
= klE~Pn
-
(9)
Ic6G~[EpPn]
Messenger gNA dM dt
= k6GE[EvPn] -
kTBM
-
k12M + ]qoN[CPa]
(10)
Template dN d---[ = k 7 B M
-
kloN[CPa]
-
(11)
k28N
Amino Acid Pool dP a d t = lc2°P*E -
(12)
kgCP~
Transport gNA-Amino Acid Complex d[ePa] dt
= kgCP~ -
kl°N[CP~]
-
(13)
k~INp[CPa]
Total Protein dE d---[ = k l ° N [ C P ~ ]
-
k16E -
k22EP~ + k 3 2 [ E P i ] -
k2EB'
(14)
Protein-Intercellular Metabolic Complex d[EP,] dt
RNA
= k22EP* -
(15)
k32[EPi]
Polymerase Cycle
Messenger RNA dMp dt
= ksGpP~ -
ksBMp
+ kllNp[CP~]
-
k13Mp
(16)
S I M U L A T I O N O F C E L L BEI-IAVIOI~
433
Template dNv dt
= ksBMp
S t a t e e q u a t i o n s o f the cell.
-
kIzNp[CP~] -
k2gNp
(1.7)
F o r the sake of conciseness, the above equations
will be described by: = 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 eel1 given by: xl = Ev
x7 = M
x13 = E
x2 = P ~
xs = M v
x l ~ = P~
x3 = [E~,P,~]
x9 = .N
x l ~ = [EvC]
x~ = B '
xlo = N v
x16 = [ E v B ]
x5 = B
xll = P~
x17 = [EP~]
x6 = C
x12 = [ C P a ]
(19)
where u e R ~ are the external inputs to the cell given by: ul = Pe
u 3 = G~
u 2 = GE
u4 = G B
u 5 = Gv
(20)
and which m u s t be determined, where k e R 29 is a vector of u n k n o w n rate constants for the cell which m u s t be determined, and x(0) e 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, a n d m u s t 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 G r o w i n g 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 w h e n a cell divides are n o t well understood.). Consider the system (18); the following division hypothesis is made: Division hypothesis. L e t t* e [0, ~ ) be the time which globally minimizes the criterion function J ( t ) ~ IIx(t) - 2x(0)l I. Then if g ( t * ) < IIx(0)[I, cell division takes place at t = t* ; i f g ( t * ) = llx(0) ll, cell division does n o t t a k e place. I f cell division takes place, the initial condition of the new generated cell is equal to ½x(t*). Figure 2 illustrates this division hypothesis. The m o t i v a t i o n for such a hypothesis is t h a t a cell in steady-state behavior, i.e. w i t h identical d y n a m i c behavior existing from one generation of cell growth to another, will have the property t h a t g ( t * ) = 0 < IIx(0)ll, and x(0) = ½x(t*). L e t x~, t~*, i = 1, 2, 3 , . . . denote the state and division time respectively of
434
E . J. D A V I S O N
,~.(o,[,
~/
0
K
(a)
IIx_{tl-2~(o)11~ lI.x(o)ll~ o (b)
Figure 2.
Criterion for deciding when
mitosis takes place. (a) Cell divides at time t* in this case. (b) Cell does not divide in this case t h e p a r e n t cell, 1st d a u g h t e r cell, 2nd d a u g h t e r cell etc. T h e n t h e behavior of a growing cell is thus described b y the following sequence of differential equations: xl = f(xl, u , k ) , i;i+ z
--'~
f(x,+ 1, u, k),
xl(0) = x °,
xt+l(0 ) -- ½x,(t*),
0 < t
--
i + 1 ,
(21) i
~
1, 2, 3, •
"
"
where t~*, i = 1, 2, 3 . . . . is f o u n d so t h a t rain t, eto, ~)[Ix,(t,) - 2x~(o)ll occurs at t~ = tt*, assuming t h a t t h e constraint: -
2x,(o)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 n o t t a k e place for i >__I, a n d the cell stops dividing. Steady-state behavior in the cell is said to h a v e t a k e n place if t h e r e exists a scalar t~* > 0 such t h a t : lira t* = t* t ~ ~o
lim Ilxi+~(t) -- xi(t)l I = O,
(23)
Vt E[O, t*]
i--* co
t~* is called the steady-state division time of t h e cell in this case. I t will be assumed in the p a p e r t h a t condition (23) m u s t be satisfied for a n o r m a l growing cell or for a cell which survives the effect of an i n p u t disturbance applied to it. I t will also be assumed, for the sake o f simplicity, t h a t a cell which stops dividing, c a n n o t survive, i.e. if (22) is not satisfied for all i, t h e n the cell dies.* 2.4. Simulation of Cell Growth. I f it is assumed t h a t the p a r a m e t e r s k, u, x ° are known, t h e n the m a t h e m a t i c a l model o f a growing cell (21) is completely • N o t e t h a t t h e topic o f cell d e a t h is a v e r y c o m p l e x s u b j e c t (1V[itchison, 1965) a n d t h e f a c t t h a t a cell s t o p s d i v i d i n g does not, of course, necessarily i m p l y t h a t it is d e a d , e.g. t h e cell m a y differentiate.
SIMULATION OF CELL BEttAVIOt:¢
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 because of 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 t~K, called I~KGS) was used to integrate the equations (21). The following parameters were used in the algorithm: max step size = 0.05 see, max error = 10 -a per unit step. The solution to (21) was calculated every 0.05 sec for i = 1, 2 , . . . , 15 with t ~ [0, T], T = 100 or until max t e to, r] IIx(t) tl >- 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 t h a n a n exponential, and so the numerical solution of (18) suddenly becomes extremely large.) Typically, the equations were integrated for t ~ [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 b y using a smaller error tolerance e.g. max error = 10 .6 per unit step. (A more accurate solution could have been obtained, b u t 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 i-dimensional global minimization was carried out by an updating procedure.~ 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,
x~(t) > 0,
Vt~[0, t~*], i = 1 , 2 , 3 , . . .
(24)
where ~ > 0 means that 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 b y 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 b y (21), (22) until condition (23) was satisfied. l" U s i n g t h i s u p d a t i n g p r o c e d u r e , t h e division t i m e t* w a s 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 c a l c u l a t e d to be e q u a l to only 0.05, t h e n c o n d i t i o n (22) is n o t satisfied a n d t h e cell dies.
436
E. J . D A V I S O N
I n this identification process it was possible to t a k e a d v a n t a g e of the i n h e r e n t r e c o v e r y p r o p e r t y of the cell itself w i t h respect to the initial condition param e t e r s 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 a p p r o x i m a t e l y equal to the correct value given b y x °* say. Then, if the s y s t e m (21)-(23) really does represent a real growing cell, one can t h i n k of x ° ¢ x °* as resulting from a disturbance being applied to the cell. I n this case the cell should 'recover' after a sufficient n u m b e r o f generations h a v e elapsed, if the disturbance is n o t too severe (i.e. if x ° is 'close enough' to x°*). This means t h a t the initial conditions of the d a u g h t e r cells will e v e n t u a l l y a p p r o a c h the correct value x °* after a sufficient n u m b e r of generations h a v e elapsed. This is illustrated in the following table (Table II) where TABLE I I Table Showing How Initial Conditions of Cell States Were Identified Parent Cell
Initial Conditions Xl(0) Assumed (Listed from x I to 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
ist Daughter Cell
Initial Conditions x2(0) Obtained (Listed from x I to x17)
e
(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 Daughter Cell (t 2 = 5.25)
Initial Conditions x3(0) Obtained (Listed from x I to x17 ) 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 (t9 : 5.25)
Initial Conditions Xl0(O) 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
l.iS
1.03
0.241
0.850
0.0115
0.0198
0.333
SIMULATION O F 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/lOth of its correct value. In this case, 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. Results 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
Pn
Xl(O) Ix2(O) x3(O) 0.170. 1.36 0.384
Pi
M N Np Pa ![CPa] E P x4(O).Xs(O)x6(O) x7(O) x8(O) x9(O) iXl0 (0) Xll(0) x12(0} x13(0: 0.218 0.304 0.114 0.334 0.244 0.078 0.066 1.14 1.03 0.240
[EpPn] B'
[EpC]
C
B
, ),I
[EpB] [EPi]]
X14 (0) xi5(0) x16(0)x17(0) 0.850
0~0114 0,0197 0.332
TABLE IV Kate Constants and External Inputs Obtained for Normal Cell (ts* = 5.25) kI
k2
k3
k4
k5
k6
1.125
1.768
1,27{)
2. i00
0.847
10.93
k7
k8
k9
kl0
k12
3.178
1.053
3.408
3.833
kll 1.018
0.973
k13
k14
k16
k17
k18
0.263
0.245
1.077
i.030
0,i09
1.606
k22 0.973
k23
0.695
0.295
k28 0.017
k32
k53
GE
0.880
1.278
0.ioo
k19
k20
2.518
2.628
k29
k15
k30
0.054
0.500
GC
GB
0.100
0,100
k31 2.091 G
P 0.i00
P O
5.75
438
E.J.
DAVISOI~I
Using these parameters the following results were obtained for the simulated cell model of (21), (22) (see Table V). TABLE V Table Showing Normal Cell Reproduction P a r e n t Cell
Initial Conditions of Cell States (Listed from Xl(O ) to x17(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
0.850
0.0114
0.0197
0.332
Final Conditions of Cell States at t I = 5.25 (Listed from
x1 (tl) to XlT(ti)). 0.341
2.73
0.771
0.438
0.609
0.231
0.671
0.490
0.158
0.133
2.29
2.07
0.480
1.71
0.0231
0.0397
0.666
1st Daughter Cell
Final Conditions of Cell States at t 2 h 5.25 (Listed from x I (t2) to x17(t2) ). 0.341
2.73
0.771
0.437
0.610
0.231
0.670
0.489
0.158
0.153
2.29
2.06
0.480
1.71
0.0231
0.0397
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°/o). A plot of the transient behavior of the normal cell model for all states is given in Figures 3-9. I t should be noted that the responses obtained agree qualitatively with Heinmet's results (Heinmet, 1966), b u t differ substantially in detail. 2.7. Example~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
S I M U L A T I O N OF CELL B E H A V I O R
2.0
P~
i.5
E
1.0
0
5
I
Figure 3.
5
t(sec)
7
Variation of Pi, P~, Pn for normal cell
2.0
2
1.5
o
E
1.0
o
Figure 4.
~
~
Variation of M, M~ for normal cell
439
440
E. J. DAVISON
2.0
~.s
E
1.0
0
3
I
5
1(see)
7
Variation of E, E~ for normal cell
Figure 5.
2.0
B'
1,O
1
Figure 6.
3
5
t (sec)
7
Variation of C, B, B ' for normal cell
SIMULATION OF CELL BEHAVIOR
2.0
2
1.5
N E o
Z
1.0
o
'i
~
'"
~ t(sec)
-t '
Figure 7. Variation of Np, N for normal cell
2.0.
1.5
1.0 ¸
0
Figure 8.
t(sec)
Variation of
[CP~], [EpC] for
normal cell
441
442
E . ft. D A V I S O N
2.0
,,, 1.5 (Ep[ -8 .N_ '-6 E O Z
1.0
B)
t(sec)
F i g u r e 9.
Variation of [EP,], n o r m a l cell
[EpP,~], [EvB ] f o r
correspond to applying various disturbances to the cell (see Tables VI and VII). I n this case (Table VI) when the gene G e is reduced b y a factor of 10, cell division ceases at the 4th generation of cells and the cell eventually dies. TABLE VI Effect of Ge --> ~aGe on Initial Conditions of Succeeding Generations of Cells (Cell Dies) Division
t*
Initial
Time
Conditions
of Cell
x 1 (0)
x 2(0)
x 3(0)
x 4(0)
• ..
x16(o)
x17(o)
0.170
1.36
0.384
0.218
...
0.0197
0.332
t I = 6.60
0.146
1.92
0.482
0.312
...
0.0174
0.332
t 2 = 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.456 ......
0.00735
0.0755
t 4 = 0.05
0.069
0.418
0.109
0.229
...
0.00360
0.037{)
t 5 = 0.0S
0.037
0.209
0.0S4
0.116
....
0.000171
0.0180
-
,..
Parent
Cell
Gell Dies
-
SIMULATION OF CELL B E H A V I O R
443
TABLE VII Effect of Pe --> ~GPe on Initial Conditions of Succeeding Generations of Cells (Cell. Dies) Division
Time
Initial Conditions of Cell
xz(O)
x2(O) x3(O) x4(o)
...
Xl6(O) x17(0)
P a r e n t Cell
0.170
1.36
0.384
0.218
...
0.0197
0.332
t I = 0.05
0.090
0.680
0.185
0.110
.....
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 b y 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 above a certain minimal concentration of nutrient medium. 3. EFFECT OF Mg 2+ STARVATm~ O~ THE C E L L 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 p u t 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 xs(0) which is x~(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 x~(0) which is x~(0) = 0, and on plotting B/0.313 × 100% vs E for 0 < t < t*, t~ = 8.1, the following results were obtained (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. J. DAVISON
TABLE VIII Table Showing Steady-state Recovery of Cell Deprived of Ribosome Initial Conditions of C~lls
Cell Division Time
Xl(0)
x2(0)
x3(0)
x4(0)
x5(0)
x16 (0)
x17 (0)
Parent Cell
0.170
1.36
0.384
0.218
0
0.0197
0.332
t I = 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
t3 = 5.00
0.169
1.43
0.400
0.218
0.313
0.0200
0.352
= 5.00
0.169
1.43
0.400
0.217
0,313
0.0201
0.351
t
t4
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, coli b u t only of a representative cell.
o
%
../t
Ribosome Content
(B)
40
/,
,H
s s '~
ao-/~'2" 0.2
0.4
0,6
0.8
Rote of Protein Synthesis (arbitrary units)
Figure 10. Comparison of experimental a n d theoretical recovery of cell from Mg 2 + starvation. (a) C o m p u t e d curve o b t a i n e d from cell model. (b) E x p e r i m e n t a l curve
4. EFFEGT OF ~ONIZINfi I~£DIATION 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
S I M U L A T I O N OF CELL B E H A V I O R
445
results for the case of single-dose and split-dose ionizing radiation disturbances. I n 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 by introducing a radiation effect hypothesis.
Radiation effect hypothesis. A dose D of ionizing radiation at time t affects a cell by 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 t h a t experimental evidence has indicated t h a t ionizing radiation actually destroys part of the various components of the cell being radiated (Little, 1968). I t will be assumed then t h a t 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 t h a t 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. TABLEIX Ionization-Radiation Experiment(Single Dose Given) Radiation factor r
Division Time of Resulting Steady-state 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
0.80
7.10
0.85
6.65
0.90
6.30
0.95
5.75
1.00 normal initial conditions
5.25 normal cell
446
E. J. DAVISON
TABLE
X
Some Typical Results Obtained for Single Dose of Ionizing Radiation Experiment Radiation Factor
r - 0,2
r = 0.6
Division Time
Initial Conditions of Cell State (Listed in order Xl(O ) to XlT(O))
t[ = 0.05
iCell does not divide.
t2 = 0.05
I
t I = 14.1
0.176 0.520 0,166 0.240 0.185 0.072 0.164
Cell dies.
0.032 0.027 O.SlO 0.750 0.084 0.70
0.116
0.0081 0.0144
0.107 t~
= O.OS Cell does not divide.
Cell dies.
t 3 = 0.05
r = 0.65
t I = 14.1
0.182 0.790 0.253 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~ = 9.95 t4 = 9.85 t~ = 9,75 Steady-state Behaviour
I
Reached
to = 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 O. 193
t[ : 9.10
r = 0.75
0.180 0.920 0.285 0.229 0.253 0.098 0.259
0.189
0.058 0.049 0.835 0.905 0.160 0.805 0.0108 0.0184 0.215 Steady-state 8ehavJour Reached
)
t 2 = 8.10 t3
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.0177 0.220
I t is seen t h a t if t h e radiation dose is t o o 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) the cell s u r v i v e s 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). I t m a y be observed in Table X t h a t for r = 0.6, t h e cell divides once before it dies, whereas for r = 0.2, the cell dies
S I M U L A T I O N O F C E L L BEI-IAVIOI~
447
immediately without dividing. This observation agrees with the following remark taken from Little (1968): "A rather perplexing feature t h a t 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 ½ of the radiation dose at t = 0 and the other ½of 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 compaxed 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. i
r{ =0,75 r2= 0,75
divisiontime = 5.;' -,i
Of RadiationSPlitDose.
{"
{ {
Deod I
0
•
9'6.B---3
sec
~
6.6
.~1/" 1/~ /+/1"// , / ~I
~ / . 4~ A 6 2 ,x Time
=
/ / / / /2/ j" ~"///~
, V/.V..:;E///.V~ ,V//~+, ,o ,2 ,4 ,6 ~6 2o
B Interval Between Doses
(See)
Figure 11. 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 rl, r2. Normal steady-state division time of the cell is 5.2 sec 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 t h a t all cell parameters are normal, except for the initial conditions which axe all multiplied by the radiation factor rl = 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 t h a n 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 . ft. D A V I S O N
the two doses of radiation are very close together, and so the overall effect of radiation approximately corresponds to a r a d i a t i o n factor of 0.75 x 0.75 = 0.56, which from the single dose experiment (Table IX) is k n o w n 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 m u c h later time ~ > 2.5, and at this time i the n a t u r a l recovery (growing) properties of the cell have started to minimize the effects of the 1st dose o n t h e cell; in particular, the m a g n i t u d e of the states of the cell are m u c h 0.05 DI = 730 rods O2=730 rods division time = 5,5 hrs.
0.04
0.03
J ~ ~ ~ / , / ~ %
Fraction Surviving Split Dose 002
O.O, , ~ / f ~/ I
0
~ - - ~ j ~ ....... v ~Experimentol I
I
I
.... I
~
~
//w ,,./
predicted Result (Qualitative Only)
~
2 4 6 B I0 12 14 16 ~. Time Interval Between Doses (hrs)
IS
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, 1)2. Normal steady-state division time of the cell is 5.5 hr larger t h a n t h e y 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 d a u g h t e r 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 n o r m a l 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 y e t (9.1 < A < 10.6), the cell will eventually divide with only the 1st dose of radiation having affected it, and thus a new d a u g h t e r cell appears at time 9.1 (see Table X). I f a 2nd dose of radiation is now applied v e r y shortly after this division time (in particular if 9.1 < A < 10.6), t h e n 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 d a u g h t e r cell is in a m u c h more fragile condition t h a n the corresponding normal cell would be (in fact, the m a g n i t u d e 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 ]3EHAVIOR
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 see. Thus for A large enough (I0.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). I f a 2nd dose of radiation is now applied very shortly after this time (in particular, ff 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 conditions 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 (A > 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. Remarlc. I t 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 O~ DISTURBANCES O~ TUE CELL: 'TuMoR-CELL' 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 b y 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 :~ N o t e t h a t t h i s is n o t a l w a y s t r u e ; it d e p e n d s o n t h e relative doses rz, ra w h e t h e r t h e 2nd m i n i m u m in t h e s u r v i v a l c u r v e appears.
450
E. J. D A V I S O N
p a r a m e t e r s of t h e cell, e x c e p t for one p a r a m e t e r only, which was c h a n g e d b y a f a c t o r of 10 or ~o (see T a b l e X I ) . I n this case it can be seen t h a t for s o m e disturbances, the cell dies, for o t h e r TABLE XI S u m m a r y of Results Obtained when Initial Conditions or Parameters of Cell are Changed x I0
l)armneter Pe
*
kl0 kll k12 k13 k14 k15 k16 k17 kls
1,7
(3.9)
Steady-state Divis'£on t i ~ t s dies dies
*
dies
Gb Gc Ge G P k1 k2 k3 k4 ks k6 k7 k8 k9
xO.1
Steady-state)livision Time t s
dies dies
dies *
1.9 (1,9)
5.9 (3.6) * *
1.0 0.75
*
l.S (I.8)
dies
(i.0)
4.1 (5.1)
(0.75)
dies
dies
dies
dies 6.0 (3.7)
*
1.8 (1.8)
*
1.9 ( 1 . 9 )
*
1.1 (1.3)
dies dies
3.6 (2.1)
2.3 (1,8)
*
1.9 (2.1)
2.3 ( I . i )
* *
1.6 (i.7) 1.6 (1.6)
7.0 (2.1]
dies
2.8 (3.2)
s.s (2.3)
2.3 (2.4)
6.2 (4.6)
dies
dies
s , l (3.9) (4.o) (s.3)
2.2 (1.7)
2.8 s.2
dies
S.6 (4.3) *
1,7
(3.8)
dies
**
dies dies
k19 k20 k21 k22 k23 k2a
dies
5.6 (5.3)
5.3 (5,3)
k29 k30 k31 k32 k33
5.3 (.5.4)
5.2 (5.4)
4.4 (S.O) 7.4 (2.7)
4.4 (5.1) 5.5 (2.2) 7.4 (2.7)
5.8 (2.6)
s.s (2.6)
Xl(0) x2(O)
dies 1.2
5.3
dies S.3 (1.3)
*
12.4 (2.3)
*
1.7 (]:.9)
5.5 (1.3)
*
1.2 (1.2)
5.6 (2:2)
*
1.2 (1.2)
S.2
SIMULATION Table Ii
x3(0): x4(0) xs(0) x6(0) xT(0) x8(O)
xg(o) Xl0(0)
CELL
5,0
dies 2.7
5.2 5.0
*
1.2
*
1.2
5.5 S.l
*
1.2
5,3
*
1.6
5.2
*
1.5
2.7
5,3 5.3
2.5
5.2
1.9
5.5
dies
5.2 5.5
x12 (0) *
Xl4 (o)
1.2
BEHAVIOI~
451
(Continued)
1.9
*
Xll(0)
Xl3(o)
OF
xls(0)
*
x16(o)
*
1.2
5,2
x17(o)
*
1.9
5,2
( ) Refers to the division time achieved when the steady-state disturbed cell
has i t s perturbed parameter changed back to i t s normal value. *
Refers to the fact that cell is behaving like a fast growing tumor-cell.
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 a/so to be observed that 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 steadystate 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 xs(0) in Table ~X I I I . 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 t h a t these fast growing cells are behaving in a way which is very
452
E.J.
DAVISON
TABLE XII Summary of Results for Sensitivity Study of Some Typical Fast Growing Tumor Cells Obtained when Initial Conditions or Parameters of Cell are Changed Steady-state Division Times (t[) for Various Initial ConDisturbance Factor
ditions or Parameters Pe
Ge
kl
k2
k7
"x2(0) x6(0) x7(0) x8(0) xi5(0) x16 (0)
x I0
1.7
1.9
1.0
•0.75 i.i
1.2
1.2
1.2
1.2
1.2
1.2
xS
2.1
2.7
0.80
0.80
1.8
2.4
2.1
4*7
2.5
2.4
x4
2.5
3.0
0.85
7.*7
2.2
2.2
2.9
2.6
4.4
3.0
2.9
x5
2.8
5.6
3"0
6.5
2.6
"2.9
3.6
3.4
4.5
5.5
3.6
x.2 .
5.5
4.4
3.9
5.7
5.2
4.1
4.5
4.5
4.8
4.3
4.5
xl.S
4.0
4.8
4.4
5.2
4.0
4.8
5.0
5.0
5.0
4.9
4.9
x I.I
4.9
S.l
S.l
5.3
4.9
5.2
5.1
5.2
5.2
5.I
5.2
xl
5.3
5.5
5.3
5.3
5.3
5.5
5.3
5.3
5.3
5.3
5.5
Disturbance Factor
Division Times Gp
k21
"1.8
1.2
x0.2
3.0
2.6
x0.4
4.4
3.9
x0.1
1.9
't s ) *" Sudden change o c c u r s a t t h i s
x0.6
4.7
4.7
x0.8
5.0
S.1
x 0.9
s.2
s.2
x 0.95
5.2
5.2
xl
5.3
5.3
point.
TABLE XIII Sensitivity Study for Fast Growing Tumor Cell Resulting from Change of Initial Condition on xs(0) Disturbance Factor Steady-state
Divi-
x 10 x 9 x 8 ~ 7 x 6 . 5 x 6 . 0 x 5.5 x 5 x 4 i x 2 x 1 . 5 x 1 . 1 i x 1
i
1.2 1.2 1.3 1 : 5 1 . 5 0
sion Time (ts)
1.65
4.95 4 . 7 4.4 4 . 9
5.0
5.2
5.3
*
* Sudden change occurs at this point.
similar to tumor cell behavior. Table XV shows the way 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
r-~
o
.
.
O o
~D
.,4
L)
> hid O 00
v I
0
0
O'X
~-. " ~
~
~
x
OF
CELL
BEtIAVIOR
453
454
E.J.
DAVISON
x2
0
2
/,
I
t
x4
t
×5
o
2/
i
:~6 t
I
t
t
x8
o
I
i
I
t
0
I
t
o
[
t
2•~//D ._ xlO
xH
o
[
I
t
0
I
t
x14 1 0
I
t
xlg
0
1
t
o
i
t
I
o
t
xl? I
xl 6
0
I
t
Figure 13. Figures showing steady-state ,transient behavior of fast growing "tumor" cell resulting from disturbance x8(0)-~ xs(0) x 9. (Ordinates are n o r , realized) 5.1. Perturbations which cause Cell to Become Tumorous. The following is a s u m m a r y 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 t h a n the normal cell, (ii) the division time of the •disturbed cell is greatly speeded up, (iii) the energy requirements of the dist u r b e d cell are much greater t h a n those of the normal cell, (iv) the disturbed cell behaves in an irreversible way, i.e. it is n o t possible to apply an opposite parameter change to the disturbed cell so t h a t it returns to its normal behavior. The cell becomes 'tumor-like' when a n y one of the following changes in the normal parameters of the cell occurs: (1)
Pe-->Pe x 10
(2) (3)
Ge-->Ge × 10 Gv --> Gv × ~-6
(4)
P~(O) --> P~(O) × 10
SIMULATION OF CELL B E H A V I O R
455
(5) [EpPn](0 ) --> [EpP~](0) x 10 (6) c(o) -~ c(o) x lO (7)
3/(0) -~ M(o) x lo
(s)
My(0 ) --> My(0) x 10
(9) N(o) ~ N(o) x lO (lO) N d o ) -+ NdO) x lO (~1)
E(o)-+ E(o) x lo
(12) [E,o](o) -~ [E~C](o) x lO (13) [EvB](0) - + [EpB](0) x 10 (14) [EP,](O) ~ [EP,](O) x 1O (15) Ep "4" Pn ~
kl
[EpPn];
/~1 ~ / d l
(16)
k2 B' + E -----> B ;
(17)
Gp + P~ k~> Gv + My;
× 10
k2 --->k2 x 10 ks ~ k5 x ~o
(18) G~ + [EvPn ] ~6> G~ + E v + M; (19) B + M kT>N;
/c6 --> k 6 x 10
k T - + k 7 x 10
(20) C + P~ k">[CP~];
kg-->ko x 10
(21) N +[CP,d klO>B + C + M + E ;
klo --> klo x 10
(22) Np + [CPa] ~1~>B + M v + C + E~,; k11--> kll x 10 (23)
Pe + E k l % p , + E ;
kls-->kla x 10
(24) E p + C . k31". ~ [EP,]; (25) E + P--~ -k'-~ [EpB]; (26) Ev + B .~---e~a" /~aa
k22 --> k22 x k2a--~ k2a x ~ .
I t would be interesting to determine those e x t e r n a l disturbances in the cell which cause the above p e r t u r b a t i o n s in the cell's p a r a m e t e r s to o c c u r - - t h i s t h e n would give insight into w h a t t y p e of e x t e r n a l disturbances in a cell cause it to b e h a v e in a t u m o r o u s way. This is a f u t u r e research area.
456
E.J.
DAVISON
6. CO~CLVSIO~S 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. (e) 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. I t 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 IIx(t*) - 2x(0)ll is minimized, provided Itx(t*) - 2x(0)II < IIx(0)II, 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)II = o.
SIMULATION OF CELL BEHAVIOR
457
This hypothesis gives a model of a cell which agrees w i t h experimental evidence for the case of a cell subject to split dose ionization radiation. F u t u r e work on the p r o b l e m would involve changing this h y p o t h e s i s a n d introducing new h y p o theses which possibly h a v e a more detailed biological basis for motivation. The above s t u d y has indicated t h a t the m a t h e m a t i c a l s t u d y of cell processes is a m e t h o d which can be used to o b t a i n insight into t h e behavior of biological processes. I t has t h e significant a d v a n t a g e t h a t simulations of complex behavior can r a p i d l y be carried o u t on a digital c o m p u t e r with all interior states of the process being monitored. F o r example, the t o t a l time required to simulate a cell dividing once was only 0.7 sec of real time on an I B M 370/165 digital c o m p u t e r ; this includes the time to plot all 17 states of the cell as a function of time a n d p e r f o r m d a t a 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 l%esearch Council of Canada under Grant A4396.
APPENDIX TERMINOLOGY AND SYMBOLS
Pools Pe = P~ = P~ = Pn =
Extraeellular nutrient pool GenerM intracellular metabolic pool Amino acid pool for protein synthesis Nucleotide pool for RNA synthesis
Enzymes E = Total protein Ep = RNA polymerase for messenger RNA (M) synthesis Genes
G~ = Genes for messenger RNA (M) synthesis Gp = Genes for messenger RNA (M~) synthesis G8 = Genes for the synthesis of I~NA fraction of ribosome
Gc = Genes for transport RNA (C) synthesis
Messengers M = Messenger (~NA) for protein (E) synthesis Mp = Messenger (I~NA) for Ep synthesis B' = RNA fraction of ribosome B = Ribosome C = Transport RNA N = Ribosome and messenger complex for protein (E) synthesis (template) Np = Ribosome and messenger complex for E~ synthesis (template) Intermediate Complexes [CPa] [EpB] [EpP,] [EP,] IN,C] 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. ~ose and E. J. Davison. 1970. "Computer Simulation of Intestinal Slow-Wave Frequency Gradient." Am. 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 Hypoxie Mammalian Cells from Sub-Lethal X - R a y Damage." Nature, Lond., 202, 1190-1193. Glueek, 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. - - . 1964b. "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: P l e n u m 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 Sy~r~posium in Math and Computer Science in the Life Sciences, Houston. Little, J . B . 1968. "Cellular Effects of Ionizing Radiation." New England J. 2~ledicine, 278, 369-376. Mallette, M. F., C. P. Clagett, A. T. Phillips and R. L. MeCarl. 1971. Introductory Biochemistry. Baltimore: Williams & Williams. McCarthy, B . J . 1962, "The Effect of Magnesium Starvation in the Ribosome Content of Escherichia Coli." Biochim Biophys Acta, 55, 880-888. Medoff, G. and M. N. Swartz. 1967. " D N A - - S t r u c t u r e and Enzyn~atic 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 R h y t h m . Bull. Math. Biology, 36, 77-89. Pollard, E. 1960. "Theoretical Aspects of the Effect of Ionizing Radiation on the Baterial Cell." Am. Naturalist, XCI¥, 71-84. Sinclair, W. K. and B. A. Morton. 1964. "Survival and Recovery in X-irradiation Synchronised Chinese Hamster Cells." Cellular Radiation Biology Proe., 18th Ann. Syrup. Fund. Ca~zcer Res., Univ. Houston, Texas. Teir, H. and T. R y t 6 m a a (eds.). 1967. Control of CeUular 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. ~]Zledicine, 276, 563-570; 276, 613-619; 276, 675-680. Weinberg, R. and B. P. Zeigler. 1969. "Computer Simulation of a Living Cell: Multilevel Control Systems." Report 08228-17-T. A n n 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. gEO]~IVED 1-17-75
VOLUME
37, 1975
S I M U L A T I O N OF CELL B E H A V I O R : NOI~MAL A N D A B N O R M A L G R O W T H
•
E. J. DAVmON 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 disturbanees, in particular when (i) a Mg~+ 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 are perturbed. The cell model results obtained agree well with experimental results for the Mga + 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.
h INTRODUCTION The purpose of this p a p e r is to find the possible m e c h a n i s m which explains w h y a n o r m a l cell s u d d e n l y becomes 'tumor-like'. This will be done b y obtaining a m a t h e m a t i c a l m o d e l of a growing cell living in a culture medium, a n d t h e n s y s t e m a t i c a l l y e x a m i n i n g the b e h a v i o r of the cell w h e n various disturbances are applied to different parts of the m a t h e m a t i c a l model of the cell. The m o t i v a t i o n for a p p l y i n g such a model a p p r o a c h is t h a t m a t h e m a t i c a l s y s t e m techniques h a v e been highly successful in explaining a n d predicting the b e h a v i o r of technical-industrial processes (see e.g. H i m m e l b l a u and Bischoff, 427
428
E . ft. D A V I S O N
1968), and recently have given insight into the mechanism of various biological processes (see e.g. Mesarovic, 1968; Diamant et al., 1970; P a l m b y 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. I-Ieinmets (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 8e, 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 t r u e 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 case 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 ceils 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 ~+ 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. ~OI~MAL 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.
E p + Pn " ~ [EpPn]
2. B' + E k-~2B k~ 3. GB + Pn _~o GB + B' k4 4. GC + Pn .-:-4. Cp + ,'.Ip kS
5. Gp + Pn -~" Gp'+ H
P
6.
k6 ~1 GE + [EpPn] ~ G]] + Ep +
7.
k7 B + !,I ~ N P
i
!o.
B kI"~4P, 1
15. C kl-~S P. 1 16. E kl-~6 p. 1.
17. g
p
k 17p.
1
18. P + E kl-~8P. + E 19. P.l + E - ~ Pn + E 20. P. + g k2-~O P
5.. P
, 1
c+P#[cPa]
lO,'N + [CPa] ~
14.
B+ C+H+E
ii. Np + [CPa] k'~11B + ~,IP + C + EP
21, Ep
+
C
k31 [EpC]
k22 22. E + Pi ~32 [EPi] 23. Ep + Bkk~3 [EpB]
12. ~Ik-~ 12 P n 13. M k-~13 P p n
+ E a
24. N-~-~ P.
1
25. N - ~ P . p I 26. p. k'-~ X 1
430
E. J . D A V I S O N
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.
RNA Synlhes~s
k9
~i
/ "-- "i
1 I l
Pool PI leal~age
7i kl
///
~<
/
/
s
INTERMEDIATE COMPLEXES
r
/
"
/ Messenll~,
ti ~~,-+-ffl-t----~
Pi
o
.7!1 -I+L
1 '
\
I
/,,.+ ..-- ~
x I +'I'~ ~-
.....
.
.
; 1<' -.*'.I,=, ~Pj
.
.
.
.
~,
t£,,.
It
.
I
=
. . . .
I
/"
1
~',.-
I
. ¢ .......................
,% POOLS Pe, PI,PI~,~o
ENZYMES Ep,E,[EPI]
MESSENGERS
Np,N,C,BiB%Mp,M
GENES GOGI~,Gp,GE
Figure 1. Flow chart showing dynamics of a cell nucleus (see Table 1) G E represents a group of genes required for protein synthesis, Gp is used for the synthesis of the enzyme E , (which is the polymerase for messenger M), gene G B produces I~NA 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 G E. Enzyme synthesis results when templates interactwith transport I~bTA and the amino acid pool complex [CPa]. The total protein E is used to convert internal pool Pi into I~NA precursors and amino acids, and to convert external pool Pe into internal pool Pi. Gene G E 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 k30), but the product is not associated with the external pool. The
S I M U L A T I O N OF C E L L B E H A V I O R
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 k~, k3~ 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 M o d e l o f 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 Fr~ction of l~ibosome dB' d--T = lc3GBP~ - lc2B'E
(1)
Ribosome dB d---{ = k 2 B ' E -
k T B M - k s B M p - k l ¢ B + kloN[CPa]
+ k l l N p [ C P ~ ] - k23E~B + k33[EpB]
(2)
Transport RNA dC
d-T = 1QGcP~ -
kgCP~ + k l ° N [ C P a ] + kllNP[CPa] - k15C
_ lc21E~C + lc3~[E,C]
(3)
General Intercellular Metabolic 1)o01 d P i = k l s E P e - k 1 9 E P i - k2oEP~ + lc3e[EP~] + k14B - ]c22.EP ~ dt
(4) + k15C + k16E + k~TE~ + ]c28N + lc29Np - k3oP~
Polymerase-Transport RNA Complex dIE, C] = lc21E, C _ k81[EpC] dt
(5)
Polymerase-RNA Complex dIE,B] dt
lc23EpB - lc33[EpB ]
(6)
432
E.J.
DAVISON
gNA Polymerase dEn = knNp[CP~] dt
-
]QTEp -
+ k33[EpB] -
k 2 1 E ~ C + ]c3a[EpC] -
k23EpB
(7)
k l E p P ~ + k6GE[E~,Pn]
Nucleotide Pool dP,~ d t = ]clgEPi -
klEpP~
ksGpP n -
+ ]c12M + k 1 3 M p -
k4GcPn
(8)
k3GBP ~
Polymerase-Nucleotide Complex d[E,P~] dt
= klE~Pn
-
(9)
Ic6G~[EpPn]
Messenger gNA dM dt
= k6GE[EvPn] -
kTBM
-
k12M + ]qoN[CPa]
(10)
Template dN d---[ = k 7 B M
-
kloN[CPa]
-
(11)
k28N
Amino Acid Pool dP a d t = lc2°P*E -
(12)
kgCP~
Transport gNA-Amino Acid Complex d[ePa] dt
= kgCP~ -
kl°N[CP~]
-
(13)
k~INp[CPa]
Total Protein dE d---[ = k l ° N [ C P ~ ]
-
k16E -
k22EP~ + k 3 2 [ E P i ] -
k2EB'
(14)
Protein-Intercellular Metabolic Complex d[EP,] dt
RNA
= k22EP* -
(15)
k32[EPi]
Polymerase Cycle
Messenger RNA dMp dt
= ksGpP~ -
ksBMp
+ kllNp[CP~]
-
k13Mp
(16)
S I M U L A T I O N O F C E L L BEI-IAVIOI~
433
Template dNv dt
= ksBMp
S t a t e e q u a t i o n s o f the cell.
-
kIzNp[CP~] -
k2gNp
(1.7)
F o r the sake of conciseness, the above equations
will be described by: = 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 eel1 given by: xl = Ev
x7 = M
x13 = E
x2 = P ~
xs = M v
x l ~ = P~
x3 = [E~,P,~]
x9 = .N
x l ~ = [EvC]
x~ = B '
xlo = N v
x16 = [ E v B ]
x5 = B
xll = P~
x17 = [EP~]
x6 = C
x12 = [ C P a ]
(19)
where u e R ~ are the external inputs to the cell given by: ul = Pe
u 3 = G~
u 2 = GE
u4 = G B
u 5 = Gv
(20)
and which m u s t be determined, where k e R 29 is a vector of u n k n o w n rate constants for the cell which m u s t be determined, and x(0) e 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, a n d m u s t 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 G r o w i n g 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 w h e n a cell divides are n o t well understood.). Consider the system (18); the following division hypothesis is made: Division hypothesis. L e t t* e [0, ~ ) be the time which globally minimizes the criterion function J ( t ) ~ IIx(t) - 2x(0)l I. Then if g ( t * ) < IIx(0)[I, cell division takes place at t = t* ; i f g ( t * ) = llx(0) ll, cell division does n o t t a k e place. I f cell division takes place, the initial condition of the new generated cell is equal to ½x(t*). Figure 2 illustrates this division hypothesis. The m o t i v a t i o n for such a hypothesis is t h a t a cell in steady-state behavior, i.e. w i t h identical d y n a m i c behavior existing from one generation of cell growth to another, will have the property t h a t g ( t * ) = 0 < IIx(0)ll, and x(0) = ½x(t*). L e t x~, t~*, i = 1, 2, 3 , . . . denote the state and division time respectively of
434
E . J. D A V I S O N
,~.(o,[,
~/
0
K
(a)
IIx_{tl-2~(o)11~ lI.x(o)ll~ o (b)
Figure 2.
Criterion for deciding when
mitosis takes place. (a) Cell divides at time t* in this case. (b) Cell does not divide in this case t h e p a r e n t cell, 1st d a u g h t e r cell, 2nd d a u g h t e r cell etc. T h e n t h e behavior of a growing cell is thus described b y the following sequence of differential equations: xl = f(xl, u , k ) , i;i+ z
--'~
f(x,+ 1, u, k),
xl(0) = x °,
xt+l(0 ) -- ½x,(t*),
0 < t
--
i + 1 ,
(21) i
~
1, 2, 3, •
"
"
where t~*, i = 1, 2, 3 . . . . is f o u n d so t h a t rain t, eto, ~)[Ix,(t,) - 2x~(o)ll occurs at t~ = tt*, assuming t h a t t h e constraint: -
2x,(o)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 n o t t a k e place for i >__I, a n d the cell stops dividing. Steady-state behavior in the cell is said to h a v e t a k e n place if t h e r e exists a scalar t~* > 0 such t h a t : lira t* = t* t ~ ~o
lim Ilxi+~(t) -- xi(t)l I = O,
(23)
Vt E[O, t*]
i--* co
t~* is called the steady-state division time of t h e cell in this case. I t will be assumed in the p a p e r t h a t condition (23) m u s t be satisfied for a n o r m a l growing cell or for a cell which survives the effect of an i n p u t disturbance applied to it. I t will also be assumed, for the sake o f simplicity, t h a t a cell which stops dividing, c a n n o t survive, i.e. if (22) is not satisfied for all i, t h e n the cell dies.* 2.4. Simulation of Cell Growth. I f it is assumed t h a t the p a r a m e t e r s k, u, x ° are known, t h e n the m a t h e m a t i c a l model o f a growing cell (21) is completely • N o t e t h a t t h e topic o f cell d e a t h is a v e r y c o m p l e x s u b j e c t (1V[itchison, 1965) a n d t h e f a c t t h a t a cell s t o p s d i v i d i n g does not, of course, necessarily i m p l y t h a t it is d e a d , e.g. t h e cell m a y differentiate.
SIMULATION OF CELL BEttAVIOt:¢
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 because of 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 t~K, called I~KGS) was used to integrate the equations (21). The following parameters were used in the algorithm: max step size = 0.05 see, max error = 10 -a per unit step. The solution to (21) was calculated every 0.05 sec for i = 1, 2 , . . . , 15 with t ~ [0, T], T = 100 or until max t e to, r] IIx(t) tl >- 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 t h a n a n exponential, and so the numerical solution of (18) suddenly becomes extremely large.) Typically, the equations were integrated for t ~ [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 b y using a smaller error tolerance e.g. max error = 10 .6 per unit step. (A more accurate solution could have been obtained, b u t 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 i-dimensional global minimization was carried out by an updating procedure.~ 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,
x~(t) > 0,
Vt~[0, t~*], i = 1 , 2 , 3 , . . .
(24)
where ~ > 0 means that 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 b y 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 b y (21), (22) until condition (23) was satisfied. l" U s i n g t h i s u p d a t i n g p r o c e d u r e , t h e division t i m e t* w a s 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 c a l c u l a t e d to be e q u a l to only 0.05, t h e n c o n d i t i o n (22) is n o t satisfied a n d t h e cell dies.
436
E. J . D A V I S O N
I n this identification process it was possible to t a k e a d v a n t a g e of the i n h e r e n t r e c o v e r y p r o p e r t y of the cell itself w i t h respect to the initial condition param e t e r s 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 a p p r o x i m a t e l y equal to the correct value given b y x °* say. Then, if the s y s t e m (21)-(23) really does represent a real growing cell, one can t h i n k of x ° ¢ x °* as resulting from a disturbance being applied to the cell. I n this case the cell should 'recover' after a sufficient n u m b e r o f generations h a v e elapsed, if the disturbance is n o t too severe (i.e. if x ° is 'close enough' to x°*). This means t h a t the initial conditions of the d a u g h t e r cells will e v e n t u a l l y a p p r o a c h the correct value x °* after a sufficient n u m b e r of generations h a v e elapsed. This is illustrated in the following table (Table II) where TABLE I I Table Showing How Initial Conditions of Cell States Were Identified Parent Cell
Initial Conditions Xl(0) Assumed (Listed from x I to 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
ist Daughter Cell
Initial Conditions x2(0) Obtained (Listed from x I to x17)
e
(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 Daughter Cell (t 2 = 5.25)
Initial Conditions x3(0) Obtained (Listed from x I to x17 ) 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 (t9 : 5.25)
Initial Conditions Xl0(O) 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
l.iS
1.03
0.241
0.850
0.0115
0.0198
0.333
SIMULATION O F 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/lOth of its correct value. In this case, 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. Results 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
Pn
Xl(O) Ix2(O) x3(O) 0.170. 1.36 0.384
Pi
M N Np Pa ![CPa] E P x4(O).Xs(O)x6(O) x7(O) x8(O) x9(O) iXl0 (0) Xll(0) x12(0} x13(0: 0.218 0.304 0.114 0.334 0.244 0.078 0.066 1.14 1.03 0.240
[EpPn] B'
[EpC]
C
B
, ),I
[EpB] [EPi]]
X14 (0) xi5(0) x16(0)x17(0) 0.850
0~0114 0,0197 0.332
TABLE IV Kate Constants and External Inputs Obtained for Normal Cell (ts* = 5.25) kI
k2
k3
k4
k5
k6
1.125
1.768
1,27{)
2. i00
0.847
10.93
k7
k8
k9
kl0
k12
3.178
1.053
3.408
3.833
kll 1.018
0.973
k13
k14
k16
k17
k18
0.263
0.245
1.077
i.030
0,i09
1.606
k22 0.973
k23
0.695
0.295
k28 0.017
k32
k53
GE
0.880
1.278
0.ioo
k19
k20
2.518
2.628
k29
k15
k30
0.054
0.500
GC
GB
0.100
0,100
k31 2.091 G
P 0.i00
P O
5.75
438
E.J.
DAVISOI~I
Using these parameters the following results were obtained for the simulated cell model of (21), (22) (see Table V). TABLE V Table Showing Normal Cell Reproduction P a r e n t Cell
Initial Conditions of Cell States (Listed from Xl(O ) to x17(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
0.850
0.0114
0.0197
0.332
Final Conditions of Cell States at t I = 5.25 (Listed from
x1 (tl) to XlT(ti)). 0.341
2.73
0.771
0.438
0.609
0.231
0.671
0.490
0.158
0.133
2.29
2.07
0.480
1.71
0.0231
0.0397
0.666
1st Daughter Cell
Final Conditions of Cell States at t 2 h 5.25 (Listed from x I (t2) to x17(t2) ). 0.341
2.73
0.771
0.437
0.610
0.231
0.670
0.489
0.158
0.153
2.29
2.06
0.480
1.71
0.0231
0.0397
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°/o). A plot of the transient behavior of the normal cell model for all states is given in Figures 3-9. I t should be noted that the responses obtained agree qualitatively with Heinmet's results (Heinmet, 1966), b u t differ substantially in detail. 2.7. Example~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
S I M U L A T I O N OF CELL B E H A V I O R
2.0
P~
i.5
E
1.0
0
5
I
Figure 3.
5
t(sec)
7
Variation of Pi, P~, Pn for normal cell
2.0
2
1.5
o
E
1.0
o
Figure 4.
~
~
Variation of M, M~ for normal cell
439
440
E. J. DAVISON
2.0
~.s
E
1.0
0
3
I
5
1(see)
7
Variation of E, E~ for normal cell
Figure 5.
2.0
B'
1,O
1
Figure 6.
3
5
t (sec)
7
Variation of C, B, B ' for normal cell
SIMULATION OF CELL BEHAVIOR
2.0
2
1.5
N E o
Z
1.0
o
'i
~
'"
~ t(sec)
-t '
Figure 7. Variation of Np, N for normal cell
2.0.
1.5
1.0 ¸
0
Figure 8.
t(sec)
Variation of
[CP~], [EpC] for
normal cell
441
442
E . ft. D A V I S O N
2.0
,,, 1.5 (Ep[ -8 .N_ '-6 E O Z
1.0
B)
t(sec)
F i g u r e 9.
Variation of [EP,], n o r m a l cell
[EpP,~], [EvB ] f o r
correspond to applying various disturbances to the cell (see Tables VI and VII). I n this case (Table VI) when the gene G e is reduced b y a factor of 10, cell division ceases at the 4th generation of cells and the cell eventually dies. TABLE VI Effect of Ge --> ~aGe on Initial Conditions of Succeeding Generations of Cells (Cell Dies) Division
t*
Initial
Time
Conditions
of Cell
x 1 (0)
x 2(0)
x 3(0)
x 4(0)
• ..
x16(o)
x17(o)
0.170
1.36
0.384
0.218
...
0.0197
0.332
t I = 6.60
0.146
1.92
0.482
0.312
...
0.0174
0.332
t 2 = 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.456 ......
0.00735
0.0755
t 4 = 0.05
0.069
0.418
0.109
0.229
...
0.00360
0.037{)
t 5 = 0.0S
0.037
0.209
0.0S4
0.116
....
0.000171
0.0180
-
,..
Parent
Cell
Gell Dies
-
SIMULATION OF CELL B E H A V I O R
443
TABLE VII Effect of Pe --> ~GPe on Initial Conditions of Succeeding Generations of Cells (Cell. Dies) Division
Time
Initial Conditions of Cell
xz(O)
x2(O) x3(O) x4(o)
...
Xl6(O) x17(0)
P a r e n t Cell
0.170
1.36
0.384
0.218
...
0.0197
0.332
t I = 0.05
0.090
0.680
0.185
0.110
.....
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 b y 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 above a certain minimal concentration of nutrient medium. 3. EFFECT OF Mg 2+ STARVATm~ O~ THE C E L L 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 p u t 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 xs(0) which is x~(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 x~(0) which is x~(0) = 0, and on plotting B/0.313 × 100% vs E for 0 < t < t*, t~ = 8.1, the following results were obtained (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. J. DAVISON
TABLE VIII Table Showing Steady-state Recovery of Cell Deprived of Ribosome Initial Conditions of C~lls
Cell Division Time
Xl(0)
x2(0)
x3(0)
x4(0)
x5(0)
x16 (0)
x17 (0)
Parent Cell
0.170
1.36
0.384
0.218
0
0.0197
0.332
t I = 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
t3 = 5.00
0.169
1.43
0.400
0.218
0.313
0.0200
0.352
= 5.00
0.169
1.43
0.400
0.217
0,313
0.0201
0.351
t
t4
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, coli b u t only of a representative cell.
o
%
../t
Ribosome Content
(B)
40
/,
,H
s s '~
ao-/~'2" 0.2
0.4
0,6
0.8
Rote of Protein Synthesis (arbitrary units)
Figure 10. Comparison of experimental a n d theoretical recovery of cell from Mg 2 + starvation. (a) C o m p u t e d curve o b t a i n e d from cell model. (b) E x p e r i m e n t a l curve
4. EFFEGT OF ~ONIZINfi I~£DIATION 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
S I M U L A T I O N OF CELL B E H A V I O R
445
results for the case of single-dose and split-dose ionizing radiation disturbances. I n 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 by introducing a radiation effect hypothesis.
Radiation effect hypothesis. A dose D of ionizing radiation at time t affects a cell by 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 t h a t experimental evidence has indicated t h a t ionizing radiation actually destroys part of the various components of the cell being radiated (Little, 1968). I t will be assumed then t h a t 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 t h a t 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. TABLEIX Ionization-Radiation Experiment(Single Dose Given) Radiation factor r
Division Time of Resulting Steady-state 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
0.80
7.10
0.85
6.65
0.90
6.30
0.95
5.75
1.00 normal initial conditions
5.25 normal cell
446
E. J. DAVISON
TABLE
X
Some Typical Results Obtained for Single Dose of Ionizing Radiation Experiment Radiation Factor
r - 0,2
r = 0.6
Division Time
Initial Conditions of Cell State (Listed in order Xl(O ) to XlT(O))
t[ = 0.05
iCell does not divide.
t2 = 0.05
I
t I = 14.1
0.176 0.520 0,166 0.240 0.185 0.072 0.164
Cell dies.
0.032 0.027 O.SlO 0.750 0.084 0.70
0.116
0.0081 0.0144
0.107 t~
= O.OS Cell does not divide.
Cell dies.
t 3 = 0.05
r = 0.65
t I = 14.1
0.182 0.790 0.253 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~ = 9.95 t4 = 9.85 t~ = 9,75 Steady-state Behaviour
I
Reached
to = 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 O. 193
t[ : 9.10
r = 0.75
0.180 0.920 0.285 0.229 0.253 0.098 0.259
0.189
0.058 0.049 0.835 0.905 0.160 0.805 0.0108 0.0184 0.215 Steady-state 8ehavJour Reached
)
t 2 = 8.10 t3
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.0177 0.220
I t is seen t h a t if t h e radiation dose is t o o 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) the cell s u r v i v e s 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). I t m a y be observed in Table X t h a t for r = 0.6, t h e cell divides once before it dies, whereas for r = 0.2, the cell dies
S I M U L A T I O N O F C E L L BEI-IAVIOI~
447
immediately without dividing. This observation agrees with the following remark taken from Little (1968): "A rather perplexing feature t h a t 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 ½ of the radiation dose at t = 0 and the other ½of 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 compaxed 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. i
r{ =0,75 r2= 0,75
divisiontime = 5.;' -,i
Of RadiationSPlitDose.
{"
{ {
Deod I
0
•
9'6.B---3
sec
~
6.6
.~1/" 1/~ /+/1"// , / ~I
~ / . 4~ A 6 2 ,x Time
=
/ / / / /2/ j" ~"///~
, V/.V..:;E///.V~ ,V//~+, ,o ,2 ,4 ,6 ~6 2o
B Interval Between Doses
(See)
Figure 11. 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 rl, r2. Normal steady-state division time of the cell is 5.2 sec 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 t h a t all cell parameters are normal, except for the initial conditions which axe all multiplied by the radiation factor rl = 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 t h a n 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 . ft. D A V I S O N
the two doses of radiation are very close together, and so the overall effect of radiation approximately corresponds to a r a d i a t i o n factor of 0.75 x 0.75 = 0.56, which from the single dose experiment (Table IX) is k n o w n 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 m u c h later time ~ > 2.5, and at this time i the n a t u r a l recovery (growing) properties of the cell have started to minimize the effects of the 1st dose o n t h e cell; in particular, the m a g n i t u d e of the states of the cell are m u c h 0.05 DI = 730 rods O2=730 rods division time = 5,5 hrs.
0.04
0.03
J ~ ~ ~ / , / ~ %
Fraction Surviving Split Dose 002
O.O, , ~ / f ~/ I
0
~ - - ~ j ~ ....... v ~Experimentol I
I
I
.... I
~
~
//w ,,./
predicted Result (Qualitative Only)
~
2 4 6 B I0 12 14 16 ~. Time Interval Between Doses (hrs)
IS
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, 1)2. Normal steady-state division time of the cell is 5.5 hr larger t h a n t h e y 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 d a u g h t e r 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 n o r m a l 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 y e t (9.1 < A < 10.6), the cell will eventually divide with only the 1st dose of radiation having affected it, and thus a new d a u g h t e r cell appears at time 9.1 (see Table X). I f a 2nd dose of radiation is now applied v e r y shortly after this division time (in particular if 9.1 < A < 10.6), t h e n 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 d a u g h t e r cell is in a m u c h more fragile condition t h a n the corresponding normal cell would be (in fact, the m a g n i t u d e 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 ]3EHAVIOR
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 see. Thus for A large enough (I0.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). I f a 2nd dose of radiation is now applied very shortly after this time (in particular, ff 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 conditions 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 (A > 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. Remarlc. I t 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 O~ DISTURBANCES O~ TUE CELL: 'TuMoR-CELL' 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 b y 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 :~ N o t e t h a t t h i s is n o t a l w a y s t r u e ; it d e p e n d s o n t h e relative doses rz, ra w h e t h e r t h e 2nd m i n i m u m in t h e s u r v i v a l c u r v e appears.
450
E. J. D A V I S O N
p a r a m e t e r s of t h e cell, e x c e p t for one p a r a m e t e r only, which was c h a n g e d b y a f a c t o r of 10 or ~o (see T a b l e X I ) . I n this case it can be seen t h a t for s o m e disturbances, the cell dies, for o t h e r TABLE XI S u m m a r y of Results Obtained when Initial Conditions or Parameters of Cell are Changed x I0
l)armneter Pe
*
kl0 kll k12 k13 k14 k15 k16 k17 kls
1,7
(3.9)
Steady-state Divis'£on t i ~ t s dies dies
*
dies
Gb Gc Ge G P k1 k2 k3 k4 ks k6 k7 k8 k9
xO.1
Steady-state)livision Time t s
dies dies
dies *
1.9 (1,9)
5.9 (3.6) * *
1.0 0.75
*
l.S (I.8)
dies
(i.0)
4.1 (5.1)
(0.75)
dies
dies
dies
dies 6.0 (3.7)
*
1.8 (1.8)
*
1.9 ( 1 . 9 )
*
1.1 (1.3)
dies dies
3.6 (2.1)
2.3 (1,8)
*
1.9 (2.1)
2.3 ( I . i )
* *
1.6 (i.7) 1.6 (1.6)
7.0 (2.1]
dies
2.8 (3.2)
s.s (2.3)
2.3 (2.4)
6.2 (4.6)
dies
dies
s , l (3.9) (4.o) (s.3)
2.2 (1.7)
2.8 s.2
dies
S.6 (4.3) *
1,7
(3.8)
dies
**
dies dies
k19 k20 k21 k22 k23 k2a
dies
5.6 (5.3)
5.3 (5,3)
k29 k30 k31 k32 k33
5.3 (.5.4)
5.2 (5.4)
4.4 (S.O) 7.4 (2.7)
4.4 (5.1) 5.5 (2.2) 7.4 (2.7)
5.8 (2.6)
s.s (2.6)
Xl(0) x2(O)
dies 1.2
5.3
dies S.3 (1.3)
*
12.4 (2.3)
*
1.7 (]:.9)
5.5 (1.3)
*
1.2 (1.2)
5.6 (2:2)
*
1.2 (1.2)
S.2
SIMULATION Table Ii
x3(0): x4(0) xs(0) x6(0) xT(0) x8(O)
xg(o) Xl0(0)
CELL
5,0
dies 2.7
5.2 5.0
*
1.2
*
1.2
5.5 S.l
*
1.2
5,3
*
1.6
5.2
*
1.5
2.7
5,3 5.3
2.5
5.2
1.9
5.5
dies
5.2 5.5
x12 (0) *
Xl4 (o)
1.2
BEHAVIOI~
451
(Continued)
1.9
*
Xll(0)
Xl3(o)
OF
xls(0)
*
x16(o)
*
1.2
5,2
x17(o)
*
1.9
5,2
( ) Refers to the division time achieved when the steady-state disturbed cell
has i t s perturbed parameter changed back to i t s normal value. *
Refers to the fact that cell is behaving like a fast growing tumor-cell.
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 a/so to be observed that 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 steadystate 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 xs(0) in Table ~X I I I . 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 t h a t these fast growing cells are behaving in a way which is very
452
E.J.
DAVISON
TABLE XII Summary of Results for Sensitivity Study of Some Typical Fast Growing Tumor Cells Obtained when Initial Conditions or Parameters of Cell are Changed Steady-state Division Times (t[) for Various Initial ConDisturbance Factor
ditions or Parameters Pe
Ge
kl
k2
k7
"x2(0) x6(0) x7(0) x8(0) xi5(0) x16 (0)
x I0
1.7
1.9
1.0
•0.75 i.i
1.2
1.2
1.2
1.2
1.2
1.2
xS
2.1
2.7
0.80
0.80
1.8
2.4
2.1
4*7
2.5
2.4
x4
2.5
3.0
0.85
7.*7
2.2
2.2
2.9
2.6
4.4
3.0
2.9
x5
2.8
5.6
3"0
6.5
2.6
"2.9
3.6
3.4
4.5
5.5
3.6
x.2 .
5.5
4.4
3.9
5.7
5.2
4.1
4.5
4.5
4.8
4.3
4.5
xl.S
4.0
4.8
4.4
5.2
4.0
4.8
5.0
5.0
5.0
4.9
4.9
x I.I
4.9
S.l
S.l
5.3
4.9
5.2
5.1
5.2
5.2
5.I
5.2
xl
5.3
5.5
5.3
5.3
5.3
5.5
5.3
5.3
5.3
5.3
5.5
Disturbance Factor
Division Times Gp
k21
"1.8
1.2
x0.2
3.0
2.6
x0.4
4.4
3.9
x0.1
1.9
't s ) *" Sudden change o c c u r s a t t h i s
x0.6
4.7
4.7
x0.8
5.0
S.1
x 0.9
s.2
s.2
x 0.95
5.2
5.2
xl
5.3
5.3
point.
TABLE XIII Sensitivity Study for Fast Growing Tumor Cell Resulting from Change of Initial Condition on xs(0) Disturbance Factor Steady-state
Divi-
x 10 x 9 x 8 ~ 7 x 6 . 5 x 6 . 0 x 5.5 x 5 x 4 i x 2 x 1 . 5 x 1 . 1 i x 1
i
1.2 1.2 1.3 1 : 5 1 . 5 0
sion Time (ts)
1.65
4.95 4 . 7 4.4 4 . 9
5.0
5.2
5.3
*
* Sudden change occurs at this point.
similar to tumor cell behavior. Table XV shows the way 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
r-~
o
.
.
O o
~D
.,4
L)
> hid O 00
v I
0
0
O'X
~-. " ~
~
~
x
OF
CELL
BEtIAVIOR
453
454
E.J.
DAVISON
x2
0
2
/,
I
t
x4
t
×5
o
2/
i
:~6 t
I
t
t
x8
o
I
i
I
t
0
I
t
o
[
t
2•~//D ._ xlO
xH
o
[
I
t
0
I
t
x14 1 0
I
t
xlg
0
1
t
o
i
t
I
o
t
xl? I
xl 6
0
I
t
Figure 13. Figures showing steady-state ,transient behavior of fast growing "tumor" cell resulting from disturbance x8(0)-~ xs(0) x 9. (Ordinates are n o r , realized) 5.1. Perturbations which cause Cell to Become Tumorous. The following is a s u m m a r y 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 t h a n the normal cell, (ii) the division time of the •disturbed cell is greatly speeded up, (iii) the energy requirements of the dist u r b e d cell are much greater t h a n those of the normal cell, (iv) the disturbed cell behaves in an irreversible way, i.e. it is n o t possible to apply an opposite parameter change to the disturbed cell so t h a t it returns to its normal behavior. The cell becomes 'tumor-like' when a n y one of the following changes in the normal parameters of the cell occurs: (1)
Pe-->Pe x 10
(2) (3)
Ge-->Ge × 10 Gv --> Gv × ~-6
(4)
P~(O) --> P~(O) × 10
SIMULATION OF CELL B E H A V I O R
455
(5) [EpPn](0 ) --> [EpP~](0) x 10 (6) c(o) -~ c(o) x lO (7)
3/(0) -~ M(o) x lo
(s)
My(0 ) --> My(0) x 10
(9) N(o) ~ N(o) x lO (lO) N d o ) -+ NdO) x lO (~1)
E(o)-+ E(o) x lo
(12) [E,o](o) -~ [E~C](o) x lO (13) [EvB](0) - + [EpB](0) x 10 (14) [EP,](O) ~ [EP,](O) x 1O (15) Ep "4" Pn ~
kl
[EpPn];
/~1 ~ / d l
(16)
k2 B' + E -----> B ;
(17)
Gp + P~ k~> Gv + My;
× 10
k2 --->k2 x 10 ks ~ k5 x ~o
(18) G~ + [EvPn ] ~6> G~ + E v + M; (19) B + M kT>N;
/c6 --> k 6 x 10
k T - + k 7 x 10
(20) C + P~ k">[CP~];
kg-->ko x 10
(21) N +[CP,d klO>B + C + M + E ;
klo --> klo x 10
(22) Np + [CPa] ~1~>B + M v + C + E~,; k11--> kll x 10 (23)
Pe + E k l % p , + E ;
kls-->kla x 10
(24) E p + C . k31". ~ [EP,]; (25) E + P--~ -k'-~ [EpB]; (26) Ev + B .~---e~a" /~aa
k22 --> k22 x k2a--~ k2a x ~ .
I t would be interesting to determine those e x t e r n a l disturbances in the cell which cause the above p e r t u r b a t i o n s in the cell's p a r a m e t e r s to o c c u r - - t h i s t h e n would give insight into w h a t t y p e of e x t e r n a l disturbances in a cell cause it to b e h a v e in a t u m o r o u s way. This is a f u t u r e research area.
456
E.J.
DAVISON
6. CO~CLVSIO~S 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. (e) 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. I t 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 IIx(t*) - 2x(0)ll is minimized, provided Itx(t*) - 2x(0)II < IIx(0)II, 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)II = o.
SIMULATION OF CELL BEHAVIOR
457
This hypothesis gives a model of a cell which agrees w i t h experimental evidence for the case of a cell subject to split dose ionization radiation. F u t u r e work on the p r o b l e m would involve changing this h y p o t h e s i s a n d introducing new h y p o theses which possibly h a v e a more detailed biological basis for motivation. The above s t u d y has indicated t h a t the m a t h e m a t i c a l s t u d y of cell processes is a m e t h o d which can be used to o b t a i n insight into t h e behavior of biological processes. I t has t h e significant a d v a n t a g e t h a t simulations of complex behavior can r a p i d l y be carried o u t on a digital c o m p u t e r with all interior states of the process being monitored. F o r example, the t o t a l time required to simulate a cell dividing once was only 0.7 sec of real time on an I B M 370/165 digital c o m p u t e r ; this includes the time to plot all 17 states of the cell as a function of time a n d p e r f o r m d a t a 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 l%esearch Council of Canada under Grant A4396.
APPENDIX TERMINOLOGY AND SYMBOLS
Pools Pe = P~ = P~ = Pn =
Extraeellular nutrient pool GenerM intracellular metabolic pool Amino acid pool for protein synthesis Nucleotide pool for RNA synthesis
Enzymes E = Total protein Ep = RNA polymerase for messenger RNA (M) synthesis Genes
G~ = Genes for messenger RNA (M) synthesis Gp = Genes for messenger RNA (M~) synthesis G8 = Genes for the synthesis of I~NA fraction of ribosome
Gc = Genes for transport RNA (C) synthesis
Messengers M = Messenger (~NA) for protein (E) synthesis Mp = Messenger (I~NA) for Ep synthesis B' = RNA fraction of ribosome B = Ribosome C = Transport RNA N = Ribosome and messenger complex for protein (E) synthesis (template) Np = Ribosome and messenger complex for E~ synthesis (template) Intermediate Complexes [CPa] [EpB] [EpP,] [EP,] IN,C] 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. ~ose and E. J. Davison. 1970. "Computer Simulation of Intestinal Slow-Wave Frequency Gradient." Am. 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 Hypoxie Mammalian Cells from Sub-Lethal X - R a y Damage." Nature, Lond., 202, 1190-1193. Glueek, 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. - - . 1964b. "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: P l e n u m 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 Sy~r~posium in Math and Computer Science in the Life Sciences, Houston. Little, J . B . 1968. "Cellular Effects of Ionizing Radiation." New England J. 2~ledicine, 278, 369-376. Mallette, M. F., C. P. Clagett, A. T. Phillips and R. L. MeCarl. 1971. Introductory Biochemistry. Baltimore: Williams & Williams. McCarthy, B . J . 1962, "The Effect of Magnesium Starvation in the Ribosome Content of Escherichia Coli." Biochim Biophys Acta, 55, 880-888. Medoff, G. and M. N. Swartz. 1967. " D N A - - S t r u c t u r e and Enzyn~atic 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 R h y t h m . Bull. Math. Biology, 36, 77-89. Pollard, E. 1960. "Theoretical Aspects of the Effect of Ionizing Radiation on the Baterial Cell." Am. Naturalist, XCI¥, 71-84. Sinclair, W. K. and B. A. Morton. 1964. "Survival and Recovery in X-irradiation Synchronised Chinese Hamster Cells." Cellular Radiation Biology Proe., 18th Ann. Syrup. Fund. Ca~zcer Res., Univ. Houston, Texas. Teir, H. and T. R y t 6 m a a (eds.). 1967. Control of CeUular 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. ~]Zledicine, 276, 563-570; 276, 613-619; 276, 675-680. Weinberg, R. and B. P. Zeigler. 1969. "Computer Simulation of a Living Cell: Multilevel Control Systems." Report 08228-17-T. A n n 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. gEO]~IVED 1-17-75