BioSystems 10 (1978) 339--347

339

© Elsevier/North-Holland Scientific Publishers Ltd.

SIMULATING FLOW MICROFLUOROMETRY EXPERIMENTS W I T H T H E C O M P U T E R L A N G U A G E CELLSIM* CHARLES E. DONAGHEY, BENJAMIN DREWINKO**, BARTHEL BARLOGIE*** and ELTON STUBBLEFIELD? Industrial Engineering Department, University of Houston, Houston, Texas 77004 and **Department of Laboratory Medicine, ***Department of Developmental Therapeutics and ?Section of Cell Biology, M.D. Anderson Hospital and Tumor Institute, Houston, Texas 77030, U.S.A.

(Received February 1E;th, 1978) CELLSIM is a digital simulation language specifically designed for simulating cell kinetics models. Recently a flow microfluorometry (FMF) command has been added to the language. In the command, the modeler indicates when his simulated FMF analyses are to take place, and what cell states are to he included. He may specify different coefficients of variance and a non-linear DNA synthesis rate for S phase cells. The FMF capability of the language has proved to he a useful tool in several applications.

1. Introduction CELLSIM is a digital simulation language specifically designed for simulating cell kinetics models. Interested readers are directed to references Donaghey and Drewinko {1975) and Donaghey (1977) for a more detailed description of the language. In developing CELLSIM every effort has been made to have the l=mguage as exportable and universal as possible. The entire interpreter has been written in standard FORTRAN, and it has been implemented on a variety of computer systems. All o u t p u t is presented on a standard line printer which has been proven to be quite, satisfactory. Recently a flow microfluorometry (FMF) c o m m a n d has been added to the language. This allows validation o f cell cycle models, because actual FMF histograms may be compared with the simulated ones obtained from CELLSIM models. When a cell population is to be analyzed by FMF equipment, the cells are stained with a flourescent dye having specific stoichiometric affinity for DNA. The stained cells are then passed, single file, through a laser beam. During the transit across the beam a fluore* This work was supported by National Cancer Institute Grants CA 19995-01 and CA ]4528.

scent light flash is emitted from each cell, the intensity of which is proportional to the a m o u n t o f fluorochrome bound to the DNA. Counters are connected to the equipment which record the number of cells having specific DNA amounts (Dean and Jett, 1974). A theoretical DNA distribution histogram for exponentially growing cells obtained by FMF is shown in Fig. 1. However, due to instrumentation error, electronic noise, light scat-

0') I0" I.l.I (,) u. 0 r,r" uJ m

=E

Z

i0 =

0 20

30

40 50 60 CHANNEL NUMBER

Fig. 1. Theoretical FMF histogram.

1 70

340

Fig. 2. Typical FMF histogram. tering, etc., the resulting FMF histogram will appear typically as shown in Fig. 2. The left hand peak on the histogram represents the cells in G1 phase (minimum amount o f DNA) and the peak to the right represents cells in G2 + M (maximum amount o f DNA). The segment between the peaks represents cells in S phase.

2. Sample model Fig. 3 shows schematically a cell model and its corresponding CELLSIM program. In this model the times in each of the phases of the cell cycle are random values from log-normal distributions. The model contains an "absorbing" state (GO), which receives 20% o f the cells exiting mitosis (M). In this model at t = 30 h, the exit from mitosis is blocked, simulating the effect of colcemid administered to the cell population at that time. The block will stay in effect until the end of the simulation at t = 50 h. The FMF command in

the program shows that FMF histograms are requested at t = 29, 32, 38 and 44 h. The command specifies that the left hand peak o f the histogram is to represent cells in the G1 and Go states of population one. The cells in S are to be in the center segment, and the right hand peak will represent cells in G2 and M. The FMF histograms obtained from this model are shown in Fig. 4. The first histogram obtained at t = 29 shows the proportion of cells in each of the 3 states (phases of the cycle) to be 0.33, 0.52 and 0.15. The coefficient of variation (CV) for the histogram is 5.0. As will be shown later, the user may specify any desired CV. The first row of numbers along the abscissa gives the channel number, and the value immediately above the channel number is the number o f cells in that channel. For example, channel 14 has 0.107 × 104 cells. The remainder of the histograms in Fig. 4 shows h o w the cells accumulate in the G~ + M peak following the block. At t = 44, 81% of the cells are in the G2 + M peak.

3. Histogram development C E L L S I M stores information on the time a cell entered a state and w h e n it is scheduled to exit. Therefore, whenever the simulation prints an F M F histogram, it is possible to find the proportion of the state each cell in that state has completed. C E L L S I M does not actually follow each cell in the simulation, since the size of a cell population m a y reach values greater than l 0 s . Instead, the system keeps information on groups o f cells, and all the cells in a group enter and leave states together during the event-based simulation. The number of groups in a model m a y grow during a simulation, and C E L L S I M has algorithms which collect and reassign groups so that m e m o r y capacity will not be exceeded. The collection algorithms are designed to have little or no disturbance on models. All cells that are inthe left peak of

341

•

8

CELL TYPES TEST; STATES(I) GI,StG2,MeGO; ABSORBING STATES[I) GO; FLOW(II GI-S(ALL}, S-G2(ALL),

G2-M(ALL),

M-GO(.2), M - G I ( . 8 ) , GI-S(ALL)I TIME IN STATES(1) S:LOGN(g,4),

GI'LOGN(-5,4),

G2:LOGN{6, l ), M:LOGN{I , I ) ; REPORTS I H O U R S ; PROLIFERATION(1) BLOCK FMF:

EXIT(1) 2q, 32,

INOCULUM(1) NO

M:2;

BETWEEN 30 AND END M; 38, 44: G I ( 1 ) + G O { I ) ,

S(II, G 2 ( 1 ) + M { I | ~

I000;

PRINT*;

SIMULATE 50 HOURS; Fig. 3. A cell kinetics model with its corresponding CELLSIM program.

the histogram are those with the m i m i m u m a m o u n t o f DNA (G1 and/ or Go states o f the model). CELLSIM finds how m a n y cells are in these states and this sum is entered into the 14th word o f a vector in the m e m o r y o f the c o m p u ter . The sum o f cells with the maxim u m a m o u n t o f DNA (G2 and/ or M states) is entered into thE: 28th word o f the vector. The thirteen words between 14 and 28 are used to store th e n u m b e r of cells that are in S phase (intermediate a m o u n t o f DNA). The cells t h a t have n o t y e t completed 1/13

o f S phase are counted and this sum placed in word 15 of the vector. Cells with more than 1/13 completed but not y e t 2/13 completed are placed in word 16. This iteration continues until all thirteen words between 14 and 28 are filled with the appropriate num ber o f cells in each segment of S phase. Each word in the vector represents one channel in the FMF equipment. The selection o f fifteen "as the n u m b e r of simulated channels produces FMF histograms on the printer close to the size obtained from t h e actual FMF assays.

342

T:29 T :58

g • I,~

I

ek; U~

g

tr

II

t~

II II II II

II II II ii

II !

!

I

X II II U

45 II II II

II II I i

0 C O00C;

I

45 II II II

It I

I

45 II li II W. .1~ II II ~, II II-I(it II rl II II II II 41- 45 II II II II 4t

11 11 II i I !

C~ 0 f5 C O

II

II I

i

II !

II i

U i

II I

II t

0 C'C C C 00

c'C I-~

w4d-

II II

II

!1 i

I

I

•If- tl l-

II

45

II II II II 45 II n

II II II 45 II II II II II II

~1 II

It

II

~1 II

i

¢ - ~ 0 0 C--C~

u'~

C"

o

(:D

u~ 'u *

UuJ

U.

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ( D O 0 0

O~

•

,-~ O4 ~e~ .,t u'~ - 0 e"- c~ O" C , ,-~ P,. 0~ ,4" u" ,C m- oL O" O ~ ~ ~

LL

U.

O O O O O O O O O O O O O O O O O O O O O O O

~_ 0

45

T: ~2 T=44 (D

45

I

45

i

4545 45

.

45 45 .

.

45 45 . . .

II II I

OOOOOO

!

II II II II i I i

II II I

II II i

O C ~O

It II t

.

.

.

45 II .

45

0

H

(D ~' ,4r

.

II

.

.

.

.

*i

.

.

II II t

II II I

tl II I

II II t

II 11 i

II II

II ,~It II I i

oCO

(DOOO

"4"

. ~

. II II !

t

i

n *

CDOOO

t

'r

•

: I

45

45

rl

45

45

}1

H

~1

i,

.

i

O ~'

~

.¢'~

g

O(:D~'OOC

¢0

COCDC

COC_,C

CC

C,C,C

C2 ¢_

C C

0

•

u.

O 0 0 0 0 0 0 0 C D o o C o 0 0 0 0 0 0 0 0 0 0

•

O OC

O O 6-5000

(.2 C O

O OOO

CO

C. O O c D ~

~ txj ,,n .,t t r , 0 e,. o. O, ¢.. ~ c.4 c~ . $ u., . ~ ~ ~C. O. CD ~ ~ j c," U

U.

u.)

Fig. 4. FMF histograms obtained for model described in Fig. 3. To the left of the histograms, are printed the time and in descending order, the computer calculated proportions of cells in G 1 +Go, S and G 2 + M phases.

4. N o r m a l i z a t i o n O n c e t h e a s s i g n m e n t o f cells is c o m p l e t e d , a n o r m a l i z a t i o n process t a k e s place. T h e h i s t o g r a m s in Fig. 4 are all printed w i t h a CV o f 5.0. This is t h e " d e f a u l t " CV i n c l u d e d in t h e program, and it m a y be o v e r r i d d e n b y t h e m o d e l e r . F o r e x a m p l e , if t h e user w a n t e d

h i s t o g r a m s w i t h a CV o f 7.0, t h e F M F c o m m a n d in t h e CELLSIM program s h o w n a b o v e w o u l d be: FMF: 2 9 , 32, 38, 44: G I ( 1 ) + G O ( l ) , S(1), G 2 ( 1 ) + M ( 1 ) : CV = 7; Using t h e CV furnished by t h e m o d e l e r , or b y

343 of a random variable lying between 13.5 and 14.5 (assuming a normal distribution with expected value 14, and standard deviation o f 0.98). The probability of cells originally assigned to channel 14 falling in channel 13, would be that o f the random variable lying between 12.5 and 13.5. These probabilities are calculated for 4SD on either side of the mean. The probability of a cell going beyond 4SD's of its original assignment is essentially zero. The cells originally assigned to channel 14 are then distributed according to these probabilities. The same normalization process is repeated for cells in each of the channels. Fig. 5 shows the effect of changing the CV on

"default", the system calculates a standard deviation for each simulated channel using the equation: (CV) (CHANNEL NUMBER)

SD =

100

For channel 14 (CV = 7.0) the standard deviation would be (7) ( 1 4 ) / 1 0 0 = 0.98. Using this SD, Channel number 14 as the mean, and assuming a normal distribution, the system calculates the distribution of the cells originally assigned to channel 14. The probability of cells originally assigned to channel 14 remaining in channel 14 would be that II

CV=5 li II ~, II II

* II II II

II II

~ II II II

II I1 !1

II ii II

'~ It II

It II II

II II II

CV=7

"

~ II ~!II ii ~ II If II

~ II II tl 11 II

TI It h II II II

I1 II

~ II

II II tl

II II tl ~

II II II

II

t~

II II II

II II II

II II II

II II II

II II II

II II II

~ II

e O e j e O O j J I I J I l J J i

Q a Q O

O

-b

II II I

II !1 q

II II II

It II II

O

Q

D

O

I

CV=5

11 tl II (I II II Q

Q

!1 ~'1 !! I ii O

f

It II II I

II II II @

I! 11 tt II II II @

J

0

II ~1~ II II tl II-~ 0

@

I

O

~

CV=9 II II -)~ It il 11 II

Ii II II II II II II N •

•

•

•

-,XII

J,f-

II II II ~ II II II ~I- ~ ~ II II II II II II II II II II II II II II II II |I II II II II II II • • • • • • .I • t, • •

it II fl II

II II

II II II II II 11 II II

I( II II II-I~II II II •

•

•

tl II II II ,XJl II II It ll'

•

•

II II II II IP II II IJ II !III II II II ~, II II II U tI i~ II II U ~I li II II Ii II II II II II II II II II II II II II II II'II II 41-

o

Fig. 5. E f f e c t o f changing c o e f f i c i e n t o f variation on the shape o f D N A histograms recorded for the same cell population.

344 the DNA histogram. Each of the histograms cdrrespond to the same cell population. Increasing the CV flattens the histogram. ,

DNA 5

5. DNA rate Another factor the modeler can control is the rate at which DNA is synthesized during S phase. If a DNA synthesis rate is not specified, the program will automatically insert a " d e f a u l t " rate. The " d e f a u l t " condition is a linear relationship between length of S phase completed and a m o u n t of DNA produced. If a user assumes that DNA is produced more rapidly in early S phase, he can specify this condition in the simulation. Fig. 6 shows such a relationship with an arbitrarily selected DNA synthesis rate. This rate can be included in an FMF command by

....

I

I

It

I

I

,

I

i

5 IO S Pig. 6. Representation o f variation in D N A synthesis Pate as a f u n c t i o n o f length in + phase. 0

furnishing points from the plot as follows: FMF: 16, 32: GI(1) + GO(I), S(1), G2(1) + M(1): DNA RATE -- 0,0 .10, .30, .30, .55, .70, .82, 1,1; For each point, the first value is the fraction II It n II It II II

.If It II |I II II

DNA

II,, II li II

DNA

;/

4f II tl II II II tl I Nil

S I

If II II II II II tl II II II II II II

41. U II I1 II II II II II II II 'K. II II 4(" "If 41' I1 II ,l II 'l 41"*

I

I

I

I

I

I

~" II 'K" II '11" II II * 4(" 'l 'l 'l I

I

I

!

4+ II |1 II It II II II 11 II II It II II II -I+ fl II t II II II II II I! 11 It I1 II I1

DNA

S

I

~I" II II II II II II

I

~. II II II II II 'l

I

4(" tl tl II II * I

i

i

i

i

!

I

i

i

DNA

1

II I| II II II II tl II II II II 41. II II II II II It II II II II ~ II II II +f 41" 11 ,l II I' II 41"

1

1

1

1

1

'15 II II "If II II tl II II II 11 11 41'

4(" II II "1" II 11 11

I

I

I

I

l

l

l

l

l

l

l

I I I

,," I1

~

41" 4{. II II II II 'if II II

~I, II tl tl

+ II II II II 41, "If 41, 41. 41, II II II II II 41, ,if II II II II II II It U II II II

~ __

M I M M I I M i M U M M M M I I i l i

• II I II I| II II I| I1 11

S

I

II tl II II II II II II It II II II II II II II II It

4+ II II II }I

II II tl II I II It II II 41. i

"II t

i

i

i

i

i

i

i

i

i

i

i

i

II 11 II II II II II n II II II II

II II ii II |1 II n

~

ii

II II II II II I! II "if II II II II II II II II 41. II N II II II N II II II NMMMMMMM

i

i

i

i

II

I I i i i i

Fig. 7. Effect of changing DNA synthesis rate on the shape of FMF histograms.

i

t

M+ i

i

i

i

II II II II II II II II II II ,If

IMMM i

i

i

i

i

i

X M i

i

i

i

i

I

345 of S phase completed, and the second value is the a m o u n t of DNA produced. Fig. 7 shows FMF histograms for the same cell population, where only the DNA synthesis rate is changed. Whenever a "non~lefault" DNA synthesis rate is furnished, the program initially assigns cells to channels by linear interpolation of the given points. It first

calculates the fraction of S phase completed by a cell, and by using this fraction, and the furnished data points, it interpolates to determine the amount of DNA synthesized. The channel to which the cell will be assigned is directly proportional to its a m o u n t of DNA. After this initial assignment of all cells, the normalization process is executed.

CELL TYPES T1; STATES(I) GIE,GIL,SE,SM,SL,G2EtG2LtMtGO; ABSORBING STATES(I) GO; TIME IN STATES(I) GIE:LOGN(?,I6), SM=LOGN(5o,I.5),

GIL=LOGNI6tlO)t SE=LOGN(I.Ot.3),

SL:LOGNI6.11.1)~G2E=LOGNI4,6)t

G2L=LOGN(2,3),

M=TRIIo3t.4~.59); INOCULUMil) 1000; FLOW(1) GIE-GIL(ALLI, GIL-SEiALL),

SE-SM(ALL), SM-SLIALL), SL-G2EiALL),

G2E-G2L(ALL), G2L-MiALL}, M-GIE(ALL); PROLIFERATION|I) M=I°6; CHANGE TIME IN STATES(1), AT TIME=41~ SE:LOGN(25.0,7) t SM=LOGN(40. O0,8), SL:LOGN(12,3|; ACCELERATE(I|t AT TIME=41~SE:.O6t SM=.125~ SL:.333; SUSPEND(1) BETWEEN 65 AND 66 SE,SMvSLtG2EtG2LtM,GIE, GIL; ACCELERATE(1)~ AT TIME=66.0ItSE:25t

SM:8, SL:3;

CHANGE TIME IN STATES(1)t AT TIME =65,SE:LOGN(1.0~o3)t

SM=LOGNI5tI°5),

SL:LOGN(6..~I.I); CHANGE TIME IN STATESI1), AT TIME=65t G2E=LOGNI6t5)t G2L=LOGNi6t3); CHANGE FLOW(1), AT TIME = 65, GIE-GIL(ALL), GIL-SEiALL), SE-SNiALL). SM-SL(ALL), SL-G2EiALL)t G2E-G2LIALL), G2L-MiALL)t M-GD|.3)~ M-GIF(°7); FMF:65.St 67, 71, 75t 78, 80, 82, 86t 93:GIEII)÷GILiII+GO(I)~ S E ( I ) - S N ( I ) - S L ( I ) ~ G2EII)+G2L(1)+M(1): DATA I:200,0~

DNA RATE=O~Oe Z , . I - O e o l ,

201t1.0;

REPORTS I HOURS; NO PRINT*; NUMBER OF GROUPS(1)=220 ; GRAPH $15 (GIE(1)+GIL(II+GO(II)/TOTAL~ SIS (SEII)+SM(I)+SLIIII/TOTAL, SIS (G2E(1)+G2LII)÷M(1))/TOTAL, SIS SE(1)/TOTAL~, TOTAL ; $IRULATE 100 HOUEI$; Fig. 8. CELLSIM Model of T, cell line with 24 h thymidine block.

Zt.7-Ot.7t I t I ;

i=======o .~==-

~o ~°

~=o

~===.

~============================= :::::::::::::::::::::::::: ::::::::::::::::::::: ~=========o

o

J

~

~°

~o

~o ~o ~~o ~==° ~.============. ~ =~==================== = = ================= = ======= = = = = = °° ~============== ° ========== * ,===.

4

============="

~===~.

~=======. ~=================. :::::::::::::::::::::::: ~==================. ~=========. ~===* ~-

O0

~.° ~° .~° ~° ~°

•

~ . = = = =t= °

.=====:

~ = = = = = = = = = = = = = = = = = = = = = =s===================== == = = == ========o= = = = = = , ~=========================

11

~° ~==========o ~========. ~===

I

~=°

~==========° ~================o

~=================° ~==============° ~=========. ~====. ~° ~° ~° ~===° ~======o

• ==========° ~================°

~="

• ° ~° ~.

s= =" ~======-

1~

::::::::::::::::::::::::

I.- °

'~

II

0

•

~° ~==o

~o~ ~. o .~

~===° ~°

~=° ~========°

~

~=========================*

~

347 6. Application of CELLSIM FMF simulation to synchronized :human l y m p h o m a cells (T1

cells) Fig. 8 describes the CELLSIM model used to simulate the behavior of T1 cells, an established human l y m p h o m a cell line, blocked by excess thymidine for 24 h. This model is based on previously published experimental data (Barlogie et al., 1976; Drewinko et al., :[977}. Assumptions used in the model will be discussed in a forthcoming publication. Fig. 9 shows actual and simulated FMF histograms at various times after the block removal. The 2 smaller peaks on the left o f the experimentally obtained histograms correspond to control cells used for quality control and are not included in the model. It is evident there is close agreement between the simulated and actual histograms.

not show directly the true characteristics of a cell population. However, the modeler can very quickly and definitely tell if the assumptions he used to develop his model were correct, hy matching simulated FMF histograms with those obtained experimentally, thus allowing iterative trials. A significant amount of insight can thus be gained in the process of developing a model, running it, studying and comparing the results, modifying it, and running it again. The T1 model shown in Fig. 5 was modified and refined over 30 times before acceptable results were obtained. Each iteration of the model generated lengthy discussions among the researchers about the fundamental behavior of the cell population and provided insight about potential avenues for experimental verification. Attention was focused on areas that would never have been considered were it not for the process of developing the model with the aid of CELLSIM.

7. Discussion The use of CELLSIM to obtain simulated FMF histograms provides researchers with an opportunity to test their assumptions about the behavior of ceil populations. The CELLSIM interpreter has been entirely written in FORTRAN and thus can be implemented on a vm~ety of computer systems. Every effort in the development of the language has been spent on ease of use and implementation. ]its aim has been to provide digital simulation capabilities to cell kinetics researchers who have little or no computer experience. The language is completely openended, in t h a t new commands and procedures may be added to the language with little or no modification of previously written code. The ease by which a CELLSIM model can be modified is, perhaps, its greatest advantage over other modeling techniques. As with any simulation model., a CELLSIM program does

References 1 Barlogie, B., B. Drewinko, T. Biichner, W. Golide, J. Schumann, W.H. Hauss and E.J. Freireich, 1976, Pulse cytophotometric analysis of cell cycle perturbation with bleomycin in vitro. Cancer Res., 36, 1176--1181. 9. Dean, P.N. and J.H. Jett, 1974, Mathematical analysis of DNA distributions derived from flow microflourometry. J. Cell. Biol. 60, pp. 523--527. 3 Donaghey, C.E., 1977, "CELLSIM II User's Manual". Department of Industrial Engineering. The University of Houston, Houston, Texas, 2nd Ed., 57 pp. 4 Donaghey, C.E. and B. Drewinko, 1975, A computer simulation program for the study of cellular growth kinetics and its application to the analysis of human lymphoma cells in vitro. Comput. Biomed. Res., 8, pp. 118--128. 5 Drewinko, B., B. Bobo, P. Roper, M.A. Malahy, B. Bariogie and B. Jansson, 1977, Analysis of the growth kinetics of a human lymphoma cell line. Proc. First Annual Meeting Cell Kinetics Society: p. 22.