A model of feedback-controlled cell populations

J. Theoret. Biol. (1968) 18, 263-279

A Model of Feedback-controlled Cell Populations JURGEN KIEFER

Radiation Center, Institute of Biophysics, Giessen,Germany (Received 31 August 1967) A generalmodelof cell populationsincluding differentiation is presented. The existenceof linear feedback-relationsbetweenthe stem-celland the mature cell demescontrolling the rate of division and differentiation is postulated. Several possibilitiesfor the elimination of mature cells are discussed.The mathematicalformalismsinvolved are outlined and it is shownthat steadystatescan be attained. The systemis investigatedunder continuousirradiation and the kinetic behaviour of a simplifiedmodel is demonstratedwith the aid of analoguecomputation. 1. Introduction During the past years many attempts have been made to describe physiological cell populations by mathematical models. They were mainly concerned with the blood forming organs, especially the erythropoietical system (seee.g. Lajtha, Oliver & Gurney, 1962; Newton, 1965). On the other hand, mathematical formalisms have been developed for the general description of such systems,e.g. the von-Foerster-equation (Trucco, 1965a,b). In this paper a general schemewill be described for a multi-deme population consisting of a stem-cell deme and of mature cell demes with the assumption of linear feedback relations between them. The term “deme” will be used in the sense as defined by Trucco (1965a). Some possibilities for regulation principles will be discussed,and for special casesthe kinetic behaviour will be outlined with the aid of analogue computation. 2. The General Scheme The general scheme of the model population is presented in Fig. 1. The central part is the stem-cell deme R from which cells may be triggered into the cell cycle, thus becoming members of T (cells in cycle) or into various possible ways of differentiation. The demes of cells in the differentiation process are termed Di, where i is the number of the special way of differentiation; mature cell demes are named Mi. R, T, Di, Mi are also used as symbols for the number of cells present in these demes.The lossesof mature 263 18 T.B.

J.KIEFER

264

cells (by death or some elimination mechanism) are described by lossfunctions Li. Cells in cycle will divide after they have completed the necessary syntheses; the time they need from the triggering for division until mitosis is called ro, the cycle time. Every cell triggered into the cycle is assumed to be doubled after r. with the two daughter cells re-entering R. No assumption will be made for the fact that the cycle time is not constant, but following a certain statistical distribution. Similarly the time necessary for the differentiation processes will be called riD and assumed to be constant.

r -r----y To

FIG.

-

R

1. Principal scheme of the model. See text for explanation.

Finally the life-times of the mature cells will be named ZiM. It should be noticed, however, that this is not in every case a reasonable value, e.g. if mature cells are eliminated at random. A further important assumption will be that there are feedback-mechanisms controlling the number of cells (these will be called “deme-sizes” in the following) in the R- and the M-demes, acting in such a manner, that the proportion of cells triggered from R into T is increased if the size of R is reduced below a certain value R,, and that the proportion of cells triggered into Di is increased if Mi is reduced below Mio. The fact that during the differentiation processes cell division may occur is accounted for by factors ci, indicating the ratio between the number of cells starting differentiation and those entering the mature cell deme. If Ci is less than unity there are cell losses during the differentiation process (“ineffective differentiation”). There are principally three possibilities for a cell entering the R-deme after mitosis. It may start a new cycle or begin differentiation or remain in a state of dormancy in respect to division or differentiation (this phase is called “Go”, Lajtha et al., 1962). Every cell has to pass this point of decision at least once during one cycle; it has been termed “dechophase” (Bullough, 1964).

MODEL

OF

FEEDBACK-CONTROLLED

CELL

POPULATIONS

265

If a system as depicted in Fig. 1 is to reach a steady state, any loss of cells must be compensated by an enhanced production of new cells. This regulation can be facilitated in three ways. The first one is to increase the fraction of cells triggered per time interval into the division cycle. This is only possible, if there is a considerable amount of cells susceptible to the triggering stimulus, that is to say, cells that are neither in the cycle nor in a differentiation process. This kind of regulation can therefore only work sufficiently if the turnover rate of the R-deme is considerably less than once per cycle-time, or expressed differently, if the number of G,-cells is not too low. The second possibility to compensate for cell losses is to shift the ratio between the numbers of cells triggered for division and differentiation towards an increased fraction of dividing cells. This mechanism will be most effective, if the turn-over rate of the R-deme is high. In this paper only the first case will be discussed in more detail; examples of the second case have been given elsewhere (Kiefer, 1966a,c). It appears to be favourable if there is only one way of differentiation, as e.g. in root meristems or the intestinal crypt cells. A third mechanism has been assumed by Newton (1965), the acceleration of the cell cycle. This case will be dealt with in the Discussion section. Let us now consider the feedback relations, indicated in Fig. 1 by dashed lines. It is assumed that as the R-deme as the M-demes have certain limiting sizes R, and Mio, respectively. If the actual values are lower, a stimulus will be given to the R-deme so that the proportion of cells triggered for division or differentiation is increased. The feedback relations may be described by functions ye(R) and ri(Mi) defined by the following relations: Number of cells triggered for division per unit time: Number

YoW.R. of cells triggered for differentiation

per unit time:

Yi(Mi1.R. The functions must fulfill the conditions: y,(R) = 0 for R = R, and 12 y,(R) > 0 for R < R, y,(M,) = 0 for Mi = Mi, and 1 2 y,(M,) > 0 for M, < Mi,. It is reasonable to assume them to be monotonic within the range considered. The simplest functions to match the given conditions are linear ones: a : constant bi: constant

266

J.

KIEFER

For this case the system of Fig. 1 may be described by the differential equations (with n being the number of differentiation ways)

g(t)=2a(l-T)

.R(t-ro)-a(l-R~) .R(t) 0

0

- i$ bi( 1 - F)

. R(t)

(la)

z(t) = ,,b,(1- M’(;-TiD)).R(t-qD)-Li. LO

Equation (lb) must be given for all n ways of differentiation (n = 4 in Fig. 1). No assumption has yet been made about the nature of the loss-functions Li. Various kinds may be constructed, but in this paper we shall constrain our considerations to three possibilities : (a) Constant efhux from Mi : Li = Ai constant. (b) Efflux proportional to the actual deme size: Li = Zi.Mi Ii: constant parameter. (c) Death of mature cells after constant life-time riM: ‘R(t-zjn-TjM). I The first question to be dealt with is whether a system as described above can attain a steady state. Assuming the existence of a steady state for t + co, one can set all derivations in respect to time equal zero and neglect the delay times. The steady-state values may then be obtained by solving the resulting algebraic equations. For the three possibilities of loss functions the following expressions are obtained : Mj(t-TjD-ZjM)

lim g = l/2(1 + W,) t-+‘w 0

(W (2b)

w, =

J

The quadratic algebraic equations have, of course, two roots, but one of them may be disregarded by considering the value for Ai = 0. In this case, R/R0 is expected to equal unity. This is only obtained with the solution given above. These considerations have also been performed for all equations to follow.

MODEL

OF

FEEDBACK-CONTROLLED

CELL

Equation (2b) holds for the ith way of differentiation. a single one, (2b) reduces to:

POPULATIONS

267

In the case of only

(24

Equations (2b) and (2~): These cases are rather similar in respect of the steady-state conditions. If the mature cells have a lifetime riM, the value of li amounts to I/riM in the steady-state. Using the abbreviation l.Mio Bi = I cibiRo

the steady-state conditions are given by

Inserting (3b) in (3a) it is obtained aR.(l-R/R,)

= R, i - buBv “cl l+B,R,/R’

(4)

Equation (4) cannot be solved generally; two special cases may be investigated : (b.1) All B, = B: then limR/R, =f(l-B+(l+B)W,) (54 t-m lim MI/M,, t-r02

= 1 - p.(l+B-(1+B)W,) t 2 Cb,

t5b)

v=l

w, =

l-4B

J

i b,/a(l+B)* v=l

(b.2) Each B, $ 1 (this does not mean necessarily that the sum of all Bi is small compared to unity). In this instance it is obtained: lim R/R, = $(l + IV,) (W Iim MJM, t-+41

= i- aB,(l-

w, =

J

l-4

W,)

2 b,B, i v=l

t6b)

i

b,B,/a. I?=1

From equations (2b), (5b) and (6b), it may be seen that the steady-state

268

J.

KIEFER

values of any M-deme are not only a function of the loss from its own population, but also of that from any other deme. The system is strongly coupled. If the sum of the losses from all demes is small compared to the maximum possible gain by division, that is to say,

Vi12 Q a& or &U4

-+ a

the steady-state values become independent :

!itRIRo= l-(~~ld,lc,)laR, lim M,/M,,

t-02

= 1 - Ai/cibiRo

U’b)

lim R/R, = 1 - i b,B,/a

t-53

(84

v=l

lim MJM,,

I’m

(74

= 1 - Bi.

@b)

The statement of the dependence of all steady-state values implies that one perturbation at any point of the system will influence the system as a whole. It should be noted that the existence of a steady state is not proved by the foregoing calculations. It is shown, however, by analogue computing (see below). The next question to be discussed is whether the transition from one steady state to another is monotonic or whether oscillations are to be expected. This problem appears to be very difficult for the general non-linear system, therefore very special assumptions will be made to get some insight into the kinetic behaviour : (1) The losses Li are to be negligible. Then the steady-state values will become R, and Mio, respectively. (2) The deviations from the steady state shall be small in each deme. This means r = (1 -R/R,) $ 1 and mi = (I -M/M,,) < 1 irrespective of time t. With the abbreviations introduced and neglecting terms quadratic in r or mi, the differential equations (1) may be written as -dr/dt(t)

= 2ar(t--o)-ar(t)-

2 b,m, i=l

-dmJdt(t) = cibimi(t-q,). The initial value for t = 0 may be sL for mi and so for r. Equation (9b) may be solved by m, = q e-Dir

Pa> Pb)

(loa)

MODEL

OF

FEEDBACK-CONTROLLED

CELL

POPULATIONS

269

where pi is given by pi = cibi i+iD, For r(t) it is obtained

fi = bici/2a @“oBi fo= &o-“ilf”

(114 WI

PO = a(2 eeporo- 1).

WC)

With the limitation given by equation (lib), equation (lob) describes the fact, that with n ways of differentiation not more than IZ extreme values are to be expected, i.e. generally the steady state will not be reached monotonically and the number of maxima and minima defines a lower limit for the number of M-demes present in the system. It must be stressed, however, that this conclusion is only valid in the case of small deviations from the steady state and for negligible losses. This result may be important for the explanation of overshoot phenomena in physiological systems. If the delay times are negligible the problem reduces to that already considered in detail for linear multi-compartment systems (Denbigh, Hicks 8z Page, 1948). 3. The System Under Continuous Irradiation As an example of constant stress exerted on the theoretical system the case of continuous irradiation will be investigated in some detail. The effect of irradiation is assumed to inhibit mitosis only in lethally damaged cells, that is to say, the flux from the T-deme to the R-deme is reduced by a factor k,(d), the surviving fraction per cell cycle time, with n being the dose-rate. It may be as low as to allow complete recovery in the resting cells. In the first example it is further assumed, that the differentiation process is not radiosensitive : R(t-L-J-U

We then find for the steady-state values under irradiation: For case (a) lim R(k,)/R, = + .(I + W,) t-+m

(124

270

J.

lim Mi(k,)/Mio t-rm

KIEFER

= 1 -Ai(2ko

WW

- l)(l - Wh)/2bici .$I $’

For case (b.1) limR(k,)/R,

(134

=+.(l-B+(l+B)W,)

t-+m n

IimMi(kO)/Mi, f’rn

Wb)

= l-a(2k,-l)(l+B-(l+B)W,))/2 2 b,, II=1 -.-____ w, = 1 - 4B i b,/a(2k, - l)( 1 + B)2 J lJ=l

For case (b.2) lim R(kJRo

= 3. (1 + W,)

t+co

lim Mi(ko)/Mio f-+02 W6

(144

= 1 - a(2k, - l)(l - We)/2 i

(14b)

b,B,

u=l

l-4

= J

i

b,B,/a(2k,

- 1)

v=l

It may be seen from equations (12a) through (14b) that real steady-state values can only be achieved if the dose-rates are not too high, or in other words, k,, not too small. No steady state can be reached in any case if k, 5 0.5. Also the final values depend strongly on the loss-functions L,; if the relative losses are high (which means high turnover-rate of the system), real steady-state values are possible only with very low-dose-rates. Let us now consider the case of the differentiation process also being radiosensitive, This is taken into account by factors ki giving the surviving fractions of differentiating cells. For the steady-state values the following expressions are yielded : For case (a) lim R(k,, k,. . . k,) = $ .(I + W,) (154 t-cc

lim M,/M,,(k,,

k, . kJ = 1 - aAi(2k, - l)(l - W,)/2biciki~~l

f--cc

w, = (bl) kl = k,.

$:" 0

(15b)

= k, = k:

lim R(k,, k)jR, = $.(l--B/k+(l+B/k)W,); ,+a,

(16a)

MODEL

OF

lim M,(k,,

t-m

FEEDBACK-CONTROLLED

CELL

= 1-a(2k,-,-l)(l+B/k-(l+B/k)Ws)/2

k)/M,

ws = J

1-4B

271

POPULATIONS

i b, “=’ (16b)

k b,/ak(2k,-l)(l+B/k)2 v=l

(b2) low dose rate: lim R/R,(k,,

t-+m

lim Mi/Mi,(k,,

1-00

kl. . .k,) = l/2.(1+

k, . . .k,) = 1-aB,(2k,-l)(l-

W,)

(1%

W,)/2ki

i o=l

w, =

b,:

’ (17b)

l-4

i b+u(‘k, - 1) ll=l Two important facts may be deduced from” equations (15a) to (17b) : (1) the steady-state values are further reduced if the differentiation is radiosensitive; (2) radiation damage in one way of differentiation will affect every other one. If we assume low dose-rates and small losses, the square root terms may be expanded, then the equations are simplified to J

lim R/R,(k,, k,. . . k,) = 1 - “& 2 (2ko- l)aR, f-+oD lim Mi/Mi,(ko, k, . . . k,) = 1 :A:jbicikiR, t-cc

lim R/R,(k,, k) = 1 -B i b,/uk(2k, ll=l t-rm lim M,/M,,(k,, k) = 1 -B/k f-rrn

- 1)

lim R/R,(k,, k, . . . k,) = 1- i b,Bk’.ju(2k, - 1) f-rc.2 u=l ” lim Mi/Mio(k,, kl . . . k,) = 1 -Bilk, 1-m

(184 (18b)

(194 Wb)

(204 VW

Two conclusions to be drawn from equations (18a) to (20b) appear to be important: (1) With low dose-rates and low turnover-rates the system is uncoupled, i.e. the different ways of differentiation are damaged independently. (2) Comparing equations (7b) with (18b) and (8b) with (19b), respectively, it may be seen that the reduction in the steady-state values of the mature

272

J. KIEFER

cell demes with very low dose-rates are only differentiation process. With low dose-rates and remain unchanged if only T-cells are damaged. Only stationary states have been investigated behaviour is also of great interest. Unfortunately (1) cannot be solved conventionally because of In order to obtain some insight into the kinetic puters is necessary. In Fig. 2 the circuit diagram

+R(I) -R(I) +I Lg,,;n ;; PM -Y

due to the damage to the low turnover-rates they will so far, although the kinetic the differential equations the delay times r,, and riD. behaviour the aid of comis given for the simulation

-R(/-R) b

To

/

FIG. 2. Circuit diagram for analogue computation L=lM.

of equation (1) with n = 2 and

of equation (1) with one way of differentiation. The delay times were generated by Pade-approximations of the second order (Giloi & Lauber, 1963); the special circuits are not given in the figure. The computations were carried out on a TELEFUNKEN RA 741 in our institute. Figure 3 shows the kinetic response after the onset of continuous irradiation for the different demes, assuming Li = Ii. Mi (case (b)) and the differentiation process being unimpaired by the radiation. The parameters used are given in the legend. It appears important to notice that the number of cells triggered for division is greatly enhanced shortly after the onset of radiation and is only reduced if the system is to run out. The nature of the regulation mechanism is very clearly shown by this feature. In Fig. 4 the respective curves are given for a system with the differentiation process also being radiosensitive (k, = k,). Here it may also be noticed that the number of cells triggered for differentiation is enhanced under irradiation.

a, 0

r0

MODEL

OF

FEEDBACK-CONTROLLED

CELL

POPULATIONS

275

4. Discussion

This paper is not aimed at comparing real systems with the theoretical one described in the foregoing sections. It is only intended to give a rough scope of the formalisms involved; a comparison between experimental facts and theoretical forecast must comprise more detailed studies of the systems under consideration, e.g. the kind of loss-functions, the number of demes, the nature of the differentiation processes, etc. Because of these reasons this section will not deal with the matching of theory and experiment, but it will try to elucidate the shortcomings and some possible advantages of the presented approach. The model may be characterized by the following features (1) It is non-linear. This is necessary to introduce feedback relations and it is possible by this means to explain how the system is able to counteract acute or continuous perturbations. (2) The model incorporates delay-times. This is necessary because of physiological grounds, as all processes, particularly the cell cycle and the differentiation processes, require time to be completed. (3) The model shows equifinality in the sense of von Bertalanffy (1953), i.e. the final steady state does not depend on the initial conditions. This statement also implies that the system is able to recover from short-time perturbations: If the values of deme-sizes reached after the perturbation being terminated are taken as initial values, the system will attain the final steadystate irrespective of the perturbation. A necessary condition, however, is that the R-deme contains still a sufficient number of cells to allow restitution. In the present state it appears very difficult to specify this statement more precisely as the equations cannot be solved analytically but it can be shown by analogue computing that after a short complete suppression of mitosis and differentiation, the final steady-state is reached again (Figs 5 and 6). For a critical discussion of the method of calculating the steady state z+de infra. (4) As indicated in the previous section the model is principally able to explain overshoot effects. Because of the mathematical difficulties involved this could not be shown generally. Unfortunately the analogue computer used was too small to investigate this problem by simulation. Thus only the very special case of small perturbations and negligible losses have been dealt with; one may speculate, however, that the possibility of the occurrence of damped oscillations may exist generally. (5) The model shows the system to act as a whole, i.e. a perturbation on one deme will generally influence all other ones. An uncoupling, however, takes place if the perturbations are small.

276

J.

KIEFER

Cd)

FIG. 5. The system when mitosis is completely suppressed (k. = 0) for a short interval. Beginning and end are indicated by arrows. Differentiation unimpaired. Parameters as in Fig. 3.

MODEL

OF

FEEDBACK-CONTROLLED

CELL

l/----t

POPULATIONS

277

(a) (b)

Time

l(----T I

Time

FIG. 6. The system when differentiation is completely suppressed (k = 0). Beginning and end are indicated by arrows. Division unimpaired. Parameters as in Fig. 3.

278

J.

KIEFER

(6) The amount of continuous perturbation (e.g. continuous irradiation) which can be tolerated depends on the loss-functions Li. If the losses are high (i.e. large A, or B,)-which means high turnover rate of the undisturbed system-only small continuous stresses are susceptible. It should be noticed that there are also some obvious shortcomings of the model proposed which may be due as to its merely hypothetical character as to its actual mathematical formulation: (I) Negative values of deme-sizes are not excluded per se, although they are biologically without sense. Therefore as an additional constraint to the model it must be stated that no deme size must be less than zero. This will certainly put additional conditions for the parameters which can be chosen and must be considered in every single case of a special system to be investigated. (2) For R > R, and Mi > M, the flow into T and Di, respectively, will become negative. This might be interpreted physiologically by a slow-down of the cell cycle or dedifferentiation, respectively, but this interpretation is rather speculative, although the model includes the possibilities. Recently Sacher & Trucco (1966) published a model of self renewing cell populations which bears some resemblance to that given here. But there are some important differences: no attempt was made to incorporate delaytimes as also particularly mentioned by the authors. No feedback control was assumed. Some of the results obtained were qualitatively similar to those yielded here although the mathematical methods were rather different. One of the first attempts to describe physiological cell systems mathematically was published by Lajtha et al. (1962) and further developed in another paper (Gilbert & Lajtha, 1965). The concept of the delay-time r0 and of the dormancy phase G, was introduced by these authors. No kind of real feedback was considered and rather complex functions had to be introduced to explain overshoot effects. The model of Iversen (1965) is based on the work of Weiss & Kavanau (1957) and Bullough (1962) and mainly concerned with the regulation in the epidermis. Only one way of differentiation is considered and only one feedback loop between division and differentiation is postulated. Furthermore it is assumed that of the two daughter cells after mitosis in every case one is triggered into a new cycle, the other starting differentiation. This scheme appears too rigid to account for the regulation in more complex systems with several ways of differentiation. Newton (1965) has suggested a model for the stem-cell kinetics of erythroand granulopoiesis in which the stabilizing feedback is assumed to act on the cell-cycle duration: a sudden removal of stem-cells is to accelerate the progression of cells in cycle towards mitosis. This concept appears difficult

MODEL

OF

FEEDBACK-CONTROLLED

CELL

279

POPULATIONS

to be widely accepted as the cell cycle duration cannot be shortened below a certain time which is necessary for the syntheses to be completed before division takes place. A model of the population kinetics in root meristems was given by the present author (Kiefer, 1966~) in which a different kind of feedback control was postulated being more appropriate to the particular system. Some features of the here outlined model have been described earlier (Kiefer, 1966a,b); particularly in the paper quoted first it was tried to compare the theory with some experimental results and to give a scheme of theoretical cell populations. Finally it may be stated that the here presented model is only to be taken as a further point of discussion in the current work of simulating physiological cell populations by mathematical models. Only a rather general line was to be drawn with no special system in mind. Applications to particular organs require necessarily more knowledge of the special features which must be obtained from experimental observations. The author is very much indebted to Professor Dr. A. Schraub, director of the Institute of Biophysics, Giessen, for steady encouragement and advice. He thanks his colleague Dr. Horst Waldvogel for stimulating discussions during the course of this investigation. REFERENCES VON BERTALANFFY,

L. (1953). “Biophysik des Fliessgleichgewichts”.

Braunschweig (Vieweg),

56pages.

(1962). Biol. Rev. 37, 307. W. S. (1964). In “Cellular Control Mechanisms and Cancer”. (P. Emmelot and 0. Miihlbock, eds), p. 124. Amsterdam: Elsevier. DENBIGH, K. G., HICKS, M. & PAGE, F. M. (1948). Trans. Faraday Sot. 44,479. GILBERT, C. W. & LAJXXA, L. G. (1965). In “Cellular Radiation Biology”, p. 474. Baltimore: The Williams & Wilkins Co. GILOI, W. & LAUBER, R. (1963). “Analogrechnen”. Berlin-Gottingen-Heidelberg: Springer, 423pages. IVERSEN,0. H. (1965). Prog. Biocyhernetics, 2,76. KIEFER, J. (1966a). Helgoliinder wiss. Meeresunters. 14, 195. KIEFER, J. (19666). Sfudia biophysics, 1, s 215. KIEFER, J. (1966c). Znt. J. Radiat. Biol. 10, 379. LAJTHA, L. G., OLIVER, R. & GURNEY, C. W. (1962). Br. J. Haemat. 8, 442. NEWTON, C. (1965). Bull. math. Biophys. 27, 275. SACHER, G. A. & TRUCCO, E. (1966). Radiat. Rex. 29, 236. TRUCCO, E. (1965a). Bull. math. Biophys. 27, 285. TRUCCO, E. (1965b). Bull. math. Biophys. 27,449. WEISS, P. & KAVANAU, J. L. (1957). J. gen. Physiol. 41, 1. BULLOUGH, BULUXJGH,

T.B.

W. S.

19