J. theor. Biol. (1991) 152, 263-270

A Discrete Model for Immune Surveillance, Tumor Immunity and Cancer DEBASHISH CHOWDHURYt~, MUHAMMAD SAHIMIt§ AND DIETRICH STAUFFERtll

t Supercomputer Center HLRZ, C / O KFA Forschungszentrum, D-5170 Jiilich, Germany, ~ School of Physical Sciences, Jawaharlal Nehru University, New Delhi 110067, India, § Department of Chemical Engineering, University of Southern California, Los Angeles, CA 90089-1211, U.S.A. and [[ lnstitut fiir Theoretische Physik, Universitiit zu K61n, D-500 K61n 41, Germany (Received on 3 August 1990, Accepted in revised form on 18 December 1990) In this paper we propose a model of tumor immunity in terms of discrete automata where each automaton describes the concentration of one particular type of cell involved in immune response. In contrast to the earlier models of normal immune response, there is more than one type of cell surface antigen in this model. As a consequence, the tumor can evade destruction through humoral response by changing its identity. However, the tumor can be killed by the killer cells through cell-mediated response unless protected by a high concentration of the suppressor T cells.

1. Introduction

Vertebrates are endowed with an i m m u n e system that patrols the lymph nodes and lymph vessels, thereby protecting the host against invading antigens, e.g. bacteria and virus ( H i l d e m a n , 1984). The i m m u n e response follows only after the antigen is detected and identified as a foreign substance. [ O f course, sometimes the i m m u n e system mistakenly identifies a part o f the host as a " f o r e i g n " substance and the i m m u n e response in such circumstances often leads to an a u t o - i m m u n e disease in the host (Cohen, 1988)]. Different types o f cells identify the antigens in different manners. For example, the B cells can identify the antigens in their native form, whereas the T cells require proper presentation o f the antigens by other cells. The response offered by various different types o f cells o f the i m m u n e system is also different. The m a c r o p h a g e s respond to all types o f antigenes, whereas the response of the lymphocytes is antigen-specific; only a specific clone o f lymphocytes fight against a specific type of antigen. There are several alternative and c o m p l e m e n t a r y mechanisms of immunity. In the so-called humoral immunity the antigens are neutralized by the antibodies produced by terminal differentiation of the B cells specifically responding to the antigen, whereas in cell-mediated immunity the T cells kill the antigen. However, it is the delicate balance between the functionings :~,§,[] Present and permanent address. 263

0022- 5193/91 / 180263+ 08 $03.00/0

O 1991 Academic Press Limited

264

D. C H O W D H U R Y

ET AL.

of the helper 1"4 cells ( T , ) and suppressor Ts cells (Ts) that maintains homeostatis, i.e. a steady state. The kinetics of immune response under various different biological circumstances have been modeled over the decades through differential equations that describe the rate of change of the concentrations of the various relevant cell types interacting among themselves (Perelson, 1988; Hoffman et al., 1988). However, often the experimental data are not accurate enough to specify the relevant ranges of the parameters involved in these differential equations. As a consequence, it is not always possible to formulate a unique model for a given immune response's kinetics. Only recently, it has been realized that simpler discrete models, with a much smaller parameter space, may often be able to describe the kinetics under several different circumstances (Atlan, 1989; Atlan & Cohen, 1989; Kaufman et al., 1985; Weisbuch & Atlan, 1988; Cohen & Atlan, 1989; Pandey et al., 1989, 1990a, b; Chowdhury & Stauffer, 1990; Neumann, 1989; Chowdhury et al., 1990; Chowdhury & Chakrabarti, 1990). Moreover, such discrete models can always be written in terms of a set of differential equations by taking the continuum limit, following well-known mathematical procedures, after specifying the relevant ranges of the parameters of the model. Such discrete models for normal immune response, auto-immune response as well as that of immune response in HIV-infected hosts have already been studied in the literature (Atlan & Cohen, 1989; Kaufman et al., 1985; Weisbuch & Atlan, 1988; Cohen & Atlan, 1989; Pandey, 1989, 1990a, b; Chowdhury & Stauffer, 1990; Chowdhury et al., 1990; Neumann, 1989; Chowdhury & Chakrabarti, 1990). However, to our knowledge, no discrete model of this type for the kinetics of t u m o u r immunity and the growth of c a n c e r o u s cells has been developed so far. In this paper we introduce such a model by extending one of the models of normal immune response proposed earlier (Chowdhury et al., 1990). This paper is organized as follows. In the next section we briefly describe the relevant aspects of tumor immunity which will be used in our model. In section 3, we describe the discrete model of tumor immunity, present the dynamical equations that describe the evolution of the system, and discuss the implications of its predictions. In the last section, we conclude by discussing how the model may be extended to more complex situations. 2. The Phenomenon of Tumor Immunity

If the rate of death of the cells is smaller than the corresponding rate of production, the tumor is called malignant and the host is said to suffer from cancer. On the one hand, the cell surface antigens are expected to be detected by the lymphocytes through a lock-and-key mechanism by which the immune system normally identifies the "foreign" substances. On the other hand, the tumor cells try to evade detection and/or avoid destruction by the immune system. Thus, the survival and growth of the tumor cells depend crucially on how the tumors evade the immune surveillance mechanism of the host and escape the killer cells as well as the antibodies. Several possible mechanisms of such evasive actions have been considered in the literature (Sell, 1984). We mention here some of the important ones that we intend to

A DISCRETE MODEL FOR IMMUNE SURVEILLANCE

265

incorporate later in our model in this paper. (i) If the tumor does not possess any cell surface antigen (i.e. if its affinity to antibodies is very low), then, obviously, the host immune system cannot detect it and respond. (ii) Even if the tumor possesses a cell surface antigen which is detected by the immune system, the tumor can survive if it can change its antigen before the antibodies or the killer cells react with it. O f course, the immune system is likely to detect the new antigen on the tumor cell surface and offer an appropriate response to the new antigen. But if the tumor can keep changing its cell surface antigens with sufficient speed, it can always escape destruction by changing its identity. (iii) However, if the tumor does not change its cell surface antigen, it can survive provided the suppressor T cells dominate the immune response, thereby protecting the tumor. As we shall show below, a set of simple dynamical equations can be formulated, by incorporating the various intercellular interactions, to describe these biological processes. The survival or destruction o f the tumor in this model depends not only on the intercellular interactions incorporated in the model, but also on the initial conditions o f the host.

3. Models of Tumor Immunity Let us begin by briefly summarizing the essential aspects o f the discrete approach.

The concentration o f each of the various cell types, rather than the individual cells, involved in the immune response is represented by a discrete two-state variable. For example, the concentration of the ith type cells is denoted by Si, where Si = 0 ("false") corresponds to the low concentration, whereas S~ = 1 ( " t r u e " ) corresponds to the high concentration o f the ith type cells. A continuous variation o f the cell concentrations cannot, of course, be described in this approach, but by assigning more than two states per variable several levels of concentration of the cells can be represented. The states o f all the variables are updated simultaneously (the so-called synchronous updating). The quantities on the right-hand side o f the dynamical equations correspond to the states of the variables at time t whereas those on the left-hand side correspond to time t + 1. The time t denotes the various stages of dynamical evolution, rather than real time in the biological world. If for all i, S~(t + 1 ) = S~(t), then the corresponding state o f the system is called a fixed point o f the dynamics. Now suppose that the antigenic epitope of type l is changed to type 2 to evade detection by the immune system as soon as T , or B cells proliferate, or when antibody concentration is high. As the host system begins to respond specifically to the antigen form 2, the tumour's antigenic epitope also begins to change its identity to a new form, say 3. If this continues the tumor gets an opportunity to always grow and evade the immune system by changing its identity. In principle, the sequence o f the various forms of the antigens can be very long. But, for simplicity o f modelling we truncate the sequence just after the second form. Nevertheless, such a simplification is not unrealistic because our model incorporates only one type of lymphocyte that responds specifically to the antigen 1, and it does not make any difference whether the antigen is in the form 2 or 3 or in any other form different from the original type 1. Moreover, it would not be difficult to extend the hierarchy by including many more forms o f antigenic epitopes, if necessary.

266

D, C H O W D H U R Y

ET AL.

We have not taken into account the killer role of the macrophages in our model, because so far as the role of the macrophages in antigen presentation is concerned, it would just delay the response by one step. We also ignore the subtle differences between the antigen recognition by the B cells and that by the T cells. In other words, we omit the details of how the antigen is detected by the T cells through presentation in the presence of the major histocompatibility complex (MHC) molecules. However, we should emphasize again that, in principle, it would pose no major difficulty to incorporate such details into the model. Among various mechanisms of tumor cell destruction, cell-mediated immunity seems to be more dominant than the humoral immunity (except when the cells are in suspension, as in e.g. leukemia) (Sell, 1984). Therefore, we formulated our model in such a way that the tumor can evade humoral immunity but can be killed by killer cells through cell-mediated immunity, unless protected by a large population of suppressor cells. Let us denote the concentrations of the antibodies, the Ts cells, the TH cells, the B cells and the killer cells, respectively, by A, S, /4, B, and K. The concentration of the tumor cells without cell surface antigens is denoted by TO, whereas those of the tumor cells with antigens of type 1 and 2 are denoted, respectively, by T1 and T2. For the three variables TO, T1 and T2, a high concentration of one excludes that of the others. For example, if most of the tumor cells have surface antigens of type 1, then, the concentration of those with no surface antigen or with surface antigen of type 2 must be relatively small. We shall take into account this mutual exclusiveness of the three variables TO, T1 and T2 in formulating the dynamical equations. Under these assumptions, the time evolution of the concentrations of the various types of cells in our model is governed by the following eight equations: A=T1.B.H,

S=H+S,

H=T1.S+H,

B=(TI+B).H

K= T1.H.S TO=TO

T I + T2

T1 = T 1 . A + H + B + K

+ TO+ T2

T 2 = ( T 2 + A + H + B ) . TO+ T1.

Here the multiplication dot stands for a logical AND, the plus sign for an OR, and the bar over a symbol for its negation, with multiplication having precedence over addition; for example, H = (T1 and not S) or/4. The first four equations are identical to those in Chowdhury et al. (1990) if we take T1 as the analog of the antigen in normal immune response. The fifth equation is also similar to the corresponding equation for concentration of the killer cells in Chowdhury et al. (1990), provided we take T1 as the analog of the (auto-) antigen. The last three equations are new. Note that there are some similarities in the structure of the last three equations. The negation of T1 + T2 in the sixth equation excludes the possibility of simultaneous high values of TO and forms T1 or T2. Similarly, parts of the last two equations also take care of the mutual exclusiveness of the variables TO, T1 and T2 mentioned earlier in this section.

A DISCRETE

MODEL

FOR

IMMUNE

SURVEILLANCE

267

The above set of eight dynamical equations can easily be solved. We find that the system possesses eleven fixed points which are listed in Table 1. The first two from the left are fixed points corresponding to a harmless state, where there are no cancerous tumor cells in the steady state. In all the other nine fixed point states malignant tumor cells exist in one form or another. Let us now discuss some o f the most interesting fixed-point aspects of this model. If the concentrations o f the eight types of cells are low, then, the concentration of none of the cells increases; this is the so-called virgin state o f the immune system. Now, suppose that the concentration o f only the Ts cells, specifically responding to T1, is high. This is the second fixed point o f Table 1 and corresponds to what has been identified as the state o f low-dose paralysis. In both these states there are no cancerous tumor cells in the host. If there is a large number o f tumor cells with cell surface antigens of type 2, then the lymphocytes specific to the antigen of type 1 do not respond to it and the state does not change with time. This is the third fixed point of our model shown in Table 1. This may be called the terminal state, since if the growth o f the tumor cells cannot be, or is not, stopped, the person with this state o f immune system will die. The effect o f the killer cells on the tumors depends crucially on the initial state o f the system. If the initial concentration o f the Ts ceils is high, then the killer cells cannot proliferate and the tumor is protected against destruction. This is the state of high-dose paralysis. According to our model this is also a fixed point state of the dynamics. This means that if the suppressor cells dominate the immune response, then the t u m o r can be protected against destruction and can continue to grow. On the other hand, if the initial concentration o f the killer cells is high and that o f the suppressors is low, then the tumor cannot escape destruction unless it senses some defensive activity of the immune system through the proliferation o f B cells, Tn cells or the antibodies. Following (Pandey, 1989, 1990a, b; C h o w d h u r y & Stauffer, 1990; Neumann, 1989; C h o w d h u r y et al., 1990) we also put our model on an L x L square lattice, with L up to 8000. Each site o f the lattice carries the eight cell variables A, S, H, B, K, TO, T1 and T2, interacting in the usual manner (see Appendix). Starting from 1% concentration o f all cell types, we end up after typically 30 iterations with an oscillation o f period two. In this final "limit cycle", we find no antibodies and no T1 cells, and all sites carry S, H, B and K. A small fraction o f sites has TO cells, TABLE 1

A list of the fixed points A S H B K TO TI T2

0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 1

0 I 1 1 0 0 0 I

0 1 1 0 0 0 0 I

0 0 0 0 0 1 0 0

0 1 l 0 0 I 0 0

0 1 1 1 0 I 0 0

0 1 0 0 0 0 0 I

0 1 0 0 0 0 I 0

0 1 0 0 0 1 0 0

268

D. C H O W D H U R Y

ET AL.

and most sites carry T2 cells. One Cray-YMP processor updated one cell variable in about a nsec. The biological significance of such lattice models is, however, controversial, since several interpretations seem possible. In these calculations we used simultaneous updating, i.e. for each iteration through all eight cell types, the input variables are those of the preceding time step. While easy to define and efficient to simulate, this method may not be biologically realistic and may lead to artificial limit cycles. Thus we also used random updating for our eight variables without a lattice; thus at each time step anew we selected by eight random numbers the type of cell to be updated next, and then used the updated value already for the next variable to be updated. We found the same fixed points; only the number of initial states which lead to fixed points without cancer cells was somewhat reduced. Similar results were found previously for related models (Strlzle, 1988; Snyder & Atlan, 1990; LeCa~r, 1989). 4. Conclusions Several models of immune surveillance and tumor immunity have been proposed earlier in the literature (de Boer & Hogeweg, 1985; Hiernaux & Lefeve, 1988; Michelson, 1988). However, in contrast to these earlier continuous models, in this paper we have introduced a simple discrete model which captures many of the crucial features of this aspect of immune response. T w o o f us (D.C, and M.S.) thank the Supercomputer Center HLRZ, Jiilich for warm hospitality, and P.V. Choudary for useful comments.

REFERENCES

ATLAN, H. (t989). Bull, math. Biol. 5, 247-253. ATLAN, H. & COHEN, R. I. (eds) (1989). Theories of Immune Networks. Heidelberg: Springer Verlag. CHOWDHURV, D. & CHAKRABARTI, B. K. (1990). Physica A 167, 628-640. CHOWDHURY, D. & STAUFFER, D. (1990). J. star Phys. 59, 1019-1042. CHOWDHURY, D., STAUFFER, D. & CHOUDARY, P. V. (1990). Z theor. Biol. 145, 207-215. COHEN, I. R. (1988). Sci. Am. 258, 34-52. COHEN, I. R. & ATLAN, H. (1989). J. Autoimmun. 2, 613-625. DE BOER, R. J. & HOGEWEG, P. (1985). In: Immunology and Epidemology (Hoffmann, G. W. & Hraba, T., eds) Lecture Notes in Biomathematics, Vol. 65. Heidelberg: Springer Verlag. HILDEMANN, W. H. (1984). In: Essentials of Immunology. Amsterdam: Elsevier. HIERNAUX, J. R. & LEFEVER, R. (1988). In: Theoretical Immunology Part 2 (Perelson, A. S., ed.) p. 19.

New York: Addison-Wesley.

HOFFMAN, G. W., KION, T. A., FORSYTH, R. B., SAGA, K. G. & COOPER-WILKIS,A. (1988). In: Theoretical Immunology Part 2 (Perelson, A. S., ed.) p. 291. New York: Addison-Wesley, KAUFMAN, M., URBAIN, J. & THOMAS, R. (1985). J. theor. Biol. i!4, 527-56t. LE CA~R, G. (1989). Physica A 157, 669-687. MlCHELSON, S. (t988). In: Theoretical Immunology Part 2 (Perelson, A. S., ed.) p. 37. New York: Addison-Wesley, NEUMAN, A. U. (1989). Physica A 162, 1-19. OLIVEIRA, P. M. C. DE (1990). In: Computing Boolean Statistical Models. Singapore: World Scientific. PANDEY, R. B. (1989). J. star Phys. 54, 997-1010. PANDEY, R. B. (1990a). J. Phys. A2 3, 4321-4331. PANDEY, R. B. (1990b)..1. star Phys. 61,235-240. PERELSON, A. S. (ed.) (1988). Theoretical Immunology Parts 1 and 2. New York: Addison-Wesley.

A D I S C R E T E M O D E L FOR I M M U N E S U R V E I L L A N C E

269

SELL, S. (1984). In: Immunology, Immunopathology and Immunity. Amsterdam: Elsevier. SNYDER, S. H. & ATLAN, H. (1990). Invariance under the order of updating in automata networks. Preprint. Sr6LZLE, S. (1988). J. star. Phys. 53, 995-1004. WEJSSUCH, G. (1991). In: Complex Systems Dynamics. New York: Addison-Wesley. WEISBUCH, G. & ATLAN, H. (1988). J. Phys. A21, LI89-192.

APPENDIX Lattice Simulation Method

In case efficient simulation of large lattices becomes biologically relevant for such immunologically motivated automata, the following FORTRAN computer program may be a useful example. It treats 64 sites in parallel with a single operation acting on 64-bit computer words, via multi-spin coding; and it is vectorized. The program of Fig. A1 can easily be modified for different Boolean models with up to N = 8 cells, by exchanging the second half of loop 7. Simpler methods using one word per site and similar methods applied to other statistical models are explained elsewhere (Weisbuch, 1991; de Oliveira, 1990). The main trick is that the expression IH.OR.IS, a logical OR of two integers IH and IS, evaluates in parallel 64 pairs of OR-conditions, corresponding to the 64 bit positions. For an L*L square lattice with L = 192, each line K is contained in only three words IN(K, J), J = 1, 2, 3. The word with J = 1 contains sites 1, 4, 7 , . . . , 190, the second word gives sites 2 , 5 , 8 , . . . , 1 9 1 , and IN(K, 3) symbolizes sites 3, 6, 9 , . . . , 192. Circular shifts in loop 5 give the left neighbors of word J = 1 and the right neighbors of word J = 3 ( = L L = L/64 in general). Two buffer lines K = 0 and K = L + 1 on top and bottom of the lattice ensure periodic boundary conditions, loop 4. The old sites IO are initialized in loops 1 and 2 with a probability p for a bit to be 1, using the random number generators RANSET and RANF. Loop 6 gives the logical OR of each site and its four neighbors (new variables IN as a function of the old variables IO); biologically this assumption is supposed to represent the spread of each cell type to empty neighbour sites. In the main loop 7 we first abbreviate the various cell types; IA means integer variable for antibody, IS for suppressor, etc. Then the interactions of the new variables IN are calculated, depending on the model, and stored in the old variables IO. Finally, loop 9 calculates with the bitcounting routine POPCNT for each cell type the number o f occupied lattice sites. On different computers the random number generators and also the logical and shift operations may have different names, like IOR(IH,IS) for IH.OR.IS.

270

1

2

4

6

7

8

9 3 10

D. C H O W D H U R Y

ET AL.

PROGRAM CANCER CELL TYPES: I = A , 2 = S , 3 = H , 4 = B , 5 = K = K I L L E R , TUMOR: 6 = T 0 , 7 = T 1 , 8 = T 2 PARAMETER ( L L = 3 , LI = LL+ 1, L=64*LL, L 2 = L + 1, N =8) D I M E N S I O N IO (0: L2, 0: L1, N), IN(0: L2, 0: Lt, N), M ( N ) DATA MAX, ISEED, P/100, 1, 0.01/ PRINT *, L, MAX, / S E E D , P CALL RANSET ( I S E E D * 2 - 1) PP1 = P + 1.0 DO 1 I = 1 , N DO 1 J = 1, LL DO1 K=I,L IO(K, J, I ) = 0 DO 2 I B I T = 0 , 63 DO21=l,N DO 2 J = 1, LL DO2K=l,L 10(K, J, I ) = IO(K, ./, I). O R . S H I F T L ( I F I X ( P P I - R A N F 0 ) , IBIT) D O 3 I T = I , MAX DO4I=I,N DO 4 J = 1, LL IO(0, J, X) = IO(L, J, t) IO(L2, J, / ) = IO(I, J, I) DO5I=I,N DO5K=l,L IO(K, 0, I) = SH1FT(IO(K, LL, I), 63) IO(K, L1, I) = SHIFT(IO(K, 1, I), 1) DO61=l,N DO6J=l, LL DO6K=I,L IN(K, J, I) = IO(K, J, I ) . O R . I O ( K - 1 , L I).OR.IO(K + 1, J, I) 1 .OR.IO(K, J - I , I).OR.IO(K, J + l , I) DO 7 J = 1, LL DO 7 K = I , L IA = IN(K, J, 1) IS = IN(K, J, 2) IH = IN(K, J, 3) I B = IN(K, J, 4) I K = IN(K, J, 5) ITO = IN(K, J, 6) IT1 = IN(K, J, 7) IT2 = IN(K, J, 8) IO(K, J, 1) = IT1.AND.IB.AND.IH 10(K, J, 2 ) = IH.OR.1S IO(K, J, 3) = (IT1.AND..NOT.IS).OR.IH IO(K, J, 4) = (ITI.OR.IB).AND.IH IO(K, J, 5 ) = IT1.AND.IH.AND..NOT.IS 10(K, J, 6) = ITO.AND..NOT.(IT1.OR.IT2) 1DUM = IH.OR.IB.OR.IA 10( K, J, 7) = IT1.AND..NOT.(IDUM.OR.ITO.OR.IT2.OR.IK) IO(K, J, 8 ) = (IT2.OR.IDUM).AND..NOT.(ITO.OR.IT1) DO81=l,N M(1) = 0 DO9 I=I,N D O 9 J = l , LL DO9K=l,L M ( I ) = M ( I ) + POPCNT(IO(K, J, I)) P R I N T 10, IT, M F O R M A T ( I X , 918) END

FnG. AI. Multi-spin coding program for eight-cell dynamics on L*L square lattice