Future Generation Computer Systems 15 (1999) 333–342

A high performance simulator of the immune response M. Bernaschi a,∗ , F. Castiglione b , S. Succi a a

b

Istituto Applicazioni Calcolo-CNR, V.le Policlinico 137, 00161 Rome, Italy Zentrum f¨ur Paralleles Rechnen (Universit¨at zu K¨oln), Wyertal 80, D-50931 Koeln, Germany Accepted 14 December 1998

Abstract The application of concepts and methods of statistical mechanics to biological problems is one of the most promising frontiers of computational physics. For instance Cellular Automata (CA), i.e. fully discrete dynamical systems evolving according to Boolean laws, appear to be extremely well suited to the simulation of the immune system dynamics. A prominent example of immunological CA is represented by the Celada–Seiden automaton that has proven capable of providing several new insights into the dynamics of the immune system response. In the present paper we describe a parallel version of the Celada–Seiden automaton. Details on the parallel implementation as well as performance data on the IBM SP2 parallel c 1999 Elsevier Science B.V. All rights reserved. platform are presented and commented on. Keywords: Immune response; Cellular automata (CA); Parallel virtual machine (PVM); Memory management

1. Introduction Most attempts to simulate the behaviour of biological systems are faced with the need of manipulating and processing vast amounts of data. The simulation of the dynamic response of the immune system (IS) is a particularly challenging task in this respect. Traditionally, the immune system response is modelled via complex networks of non-linear ordinary differential equations (ODEs). This approach is certainly valuable, but suffers from a number of drawbacks, primarily a severe exposure to numerical drifts and/or instabilities intrinsic to the floating-point representation of real numbers in digital computers. Luckily, immunological research may get significant benefits from the resort to Cellular Automata (CA) models, i.e. fully discrete dynamical systems, whose fictitious dynamics can be tuned so as to ‘mimic’ the behaviour of the real system at least on a macroscopic scale. A clear advantage of CA over ODEs is intrinsic numerical stability, which is strictly related to the boolean character of the dynamic variables [1]. One of the most prominent attempts to cope with the quest for biological fidelity is the ImmSim automaton, developed by Celada and Seiden [2,3]. ImmSim still belongs to the class of Immunological Cellular Automata (ICA), but its degree of sophistication sets it apart from simpler CA in the Ising-like class [4,5]. As a consequence, ∗

Corresponding author. E-mail: [email protected]

c 0167-739X/99/$ – see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 7 3 9 X ( 9 8 ) 0 0 0 7 8 - 8

334

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

ImmSim is very demanding from the computational viewpoint. This is the stage where parallel processing can offer the chance of a real breakthrough to immunological research. Armed with these considerations, we set out to develop a parallel version of the ImmSim automaton (ParImm). ParImm has been designed according to criteria of openness, and modularity, so as to permit smooth upgrades and addition of new features (cells, molecules, interactions, and so on) for future investigations. The paper is organized as follows: in Section 2 we give some brief explanation of the main fuctionalities of ImmSim. Section 2.2 deals with the optimization of the memory management to reduce the memory requirement and in Section 3 we discuss details of the two phase-version of the parallel code (ParImm versions 1.0 and 2.0). Finally, in Section 3.3 some data pertaining to parallel performance and results are presented and commented on. Section 4 concludes the paper with future perspectives.

2. From biology to computer modeling The IS comprises different classes of cells all cooperating to a single goal: the defence of the host organism from attacks of potentially offending invaders. In the model of the IS we use, the entities are free to move in a CA-like body rectangle, and all the interactions among them take place on-site, locally, each independently of the others. The dynamics of the system is mainly governed by these local interactions. Nothing changes in the system (apart from normal fluctuations in the population dynamics) until an antigen (Ag) (which identifies a generic substance, molecule, virus, bacteria and the like, potentially dangerous for the host organism) gets inside the body (infection). When this happens the system starts the recognition process and the response terminates when the Ag is eliminated by either humoral (antibodies (Abs)) and/or cellular response. At the same time, memory cells are produced by clonal division of lymphocites; the memory of what happened sets the system in such a way that it responds promptly to future infections by the same Ag. An infection is simulated by spreading randomly the Ag over the body rectangle. Among different kind of cells represented (B, T, AP, PLB), the real chief characters are the B cells (lymphocyte B). Upon stimulation, theses produce a clone of PLB (Plasma B cells) which will subsequently secrete Abs in charge of neutralizing the Ags. Another peculiar class of cells is the Ag processing cell (APC) or macrophage representing the non-specificity of the sytem, that is the ability to capture foreign entities without any recognition process. Actually, the APC’s just ‘pre-process’ the Ags, by splitting them in peptides which are bound to the molecule of the major histocompatibility complex (MHC). In ImmSim both cell receptors and molecules are binary bit strings of length `. The binding between two strings occurs with a probability which is a function of the matching (Hamming distance in {0, 1}` ) of the two strings. So, when cells and Ags come close together, they may interact under one condition: their receptors must have a matching value at least equal to mc (a threshold value usually called MinMatch). Every binding event induces a cell-state-change to represent the various phases of the immune response. 2.1. Model description It is important to characterize the phases of the immune system response as they are triggered by events and external conditions. The starting point is the first injection of Ags at time step t0 . This may take place anytime after the simulation starts; the system is built to maintain a steady state of the global population of cells if no infection is applied. Moreover, initially there are neither plasma cells and Abs, nor Ab–Ag complexes in the body. The various steps of the dynamic evolution are: (i) t0 : injection of Ags (the host has been infected); (ii) t1 = t0 + δ1 : B cells bind Ags with a certain affinity; AP cells bind Ags aspecifically; (iii) t2 = t1 + δ2 : B cells and AP cells process Ag internally (MHC bind Ag-peptides) ⇒ they expose the MHC/peptide groove on the surface;

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

335

Fig. 1. B and T cells.

Fig. 2. Antigens and immuno complexes.

(iv) t3 = t2 + δ3 : T cells bind AP and/or B cells which expose the MHC-complex; ⇒ both B and T cells get stimulated; (v) t4 = t3 + δ4 : stimulated T cells divide in T virgin and T memory; stimulated B cells divide in B memory and PLB cells; (vi) t5 = t4 + δ5 : PLBs secrete Abs; (vii) t6 = t5 + δ6 : Abs bind Ags. A further contact with the same Ag causes the aforementioned process to be repeated again. In Fig. 1 we report the populations of the ten receptor strings (for B, T cells) which show the best matching with the Ag. The binary strings are expressed in decimal notation in the figure-titles. Fig. 2 shows Ags and Immuno complexes (ICs), respectively. With the so called standard parameter set [2,3], the model exhibits the primary and secondary immune response as shown in Figs. 1 and 2. Here the size of the automaton is 16 × 16 and the bit string length is equal to 8. The Ags are injected at time 0 and 120 (Fig. 2). After some time (10 time steps) the clonal growth of B cells and also T cells begins to occur (Fig. 1). In Fig. 1 we see that the highest growth is for the string-population of B cells 20410 because of the affinity with the Ag-peptide 5110 = 001100112 (20410 = 110011002 is the only perfect match string of 5110 ). The other string-populations are most likely in the 7-bit match string set (205, 206, 200, 196, 220, 236, 140 and 76). Their

336

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

growth is strictly depending on the T cells growth. During the response new T and B memory and virgin cells are produced. The highest population of T cells depends on the MHC molecules. In fact they get stimulated if their receptor binds the MHC/peptide groove on the B or AP cells. String 9410 (010111102 ) is the perfect match string of 101000012 which is composed by half MHC molecule 16710 = 101001112 and half Ag-peptide 1710 = 000100012 (refer [2] for details). After the Ags disappear (about time step 30), the number of cells decreases. As to the memory cells their reduction is slower because of the longer lifetime. It is noteworthy that the time required to eliminate the Ag after the second injection is about three times shorter. 2.2. Dynamic memory allocation In a standard CA there is a simple correspondence between entities of the automaton and sites of the grid. Usually each site contains either one or zero instances of an entity and the number of possible entities is very small (in the well-known game of life there is just one possible entity). This feature makes the computation of the interactions among entities very simple. Unfortunately it imposes a severe restriction on the number of instances of each entity that can be represented (this number cannot be greater than the number of lattice points). For a biological model, this is a too strong limitation. As a matter of fact, the IS is a dynamic system with a population that, due to the combination of multiple mechanisms like evolution, learning and memory, may change significantly during the simulation. For such reason, we ended up with a computationally harder but much more flexible model. The most serious issue we had to address was the huge amount of memory required to store the information which describe the attributes of the entities. In the early stage of the translation from the original APL code into a C code we adopted a static data structure for the entities of the simulation. We used a number of multidimensional arrays, as many as the classes of entities represented in the code. The data structures were not homogeneous because of the different number of attributes of each entity. The attempt to store all the combination of attributes in a single multidimensional array by lumping together all indistinguishable cells, had good and bad points. By lumping together the cells, it is possible, to have in site X, say 1 000 517 B cells with receptor (001100112 ) all in a ‘virgin’ state not exposing Ag-peptide. The cost of such information is independent of the number of cells and requires just four bytes (i.e. the size of an integer on most machines). On the other hand we cannot assert a priori which information we need at run time because of the stochastic nature of the system. As a consequence four bytes are required for each possibility. Let’s take for example the B class of cells; we specify the following attributes: – cell position within the body ∈ {0, D − 1} – receptor ∈ {0, 2` − 1} – state ∈ {0, N States} – duplication step ∈ {0, DupSteps} – MHC/peptide grove ∈ {0, 2` − 1} this means 2` × D × N States × DupStep × 2` × sizeof(int) bytes. With this scheme for a standard 8-bit system, with fixed small Ag and just one MHC molecule, we need (` = 8, D = 16 × 15, N states = 5, DupStep = 4), 4.68 MB which is a resonable amount of memory. However, for a 12-bit system with D = 30 × 30, N states = 5, DupStep = 4 and variable (but still small) number of Agepitopes/peptides and MHC molecules, the requirement increases to 1.17 GB. It is clear that such a static structure is not affordable and that a more flexible, although more complex, dynamic memory allocation technique is in order. Having chosen the attributes for each entity the information required for each cell is organized in blocks of variables, one block for each cell. In such way each cell gets a clear identity (record) in memory and can be easily told apart from the others. The blocks are linked in lists, one for each class of cells (B list, T list, PLB list, APC list). The lists are initialized at startup time and are dynamically managed at run time. When a cell dies out it is removed from the list. Similarly when a virgin cell comes in (from the bone marrow, lymph nodes and/or for clone division of stimulated cells) memory is allocated, the variables-attributes are filled with appropriate values and the

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

337

Fig. 3. Dynamic memory requirement.

whole block becomes the head of the list. The memory grows linearly with the number of cells, following the clonal growth of cells during the immune response. As a result there are bursts in memory requirements right after the injection of antigens. In Fig. 3 it is shown the amount of memory (in bytes) required both initially and after the first and the second response as a function of the initial number of cells. 2.2.1. List permutation In the IS each cell has thousands of receptors on its surface. So it is able to bind more than one entity. Moreover, there are situations in which two cells bind a third one at the same time. With our simulator this is not possible. In the CS model, each cell has a single receptor, so the binding events follow a greedy paradigm: a successful binding event removes the tied entities from the eliciting set to prevent further interactions in the same time step. It is quite clear that such mechanism may introduce an artifical bias in the simulation. Indeed, if the scanning order during the interaction procedures were simply the cells’ lists order, the cells close to the head of the list would get more chances to catch the Ags. To alleviate such problem the list of the cells is reordered at each time step. Basically this is a two-step procedure: (i) pointers to the list elements are copied to an array of pointers in random order; (ii) a new list of cells is built with the order defined by the array of pointers. 2.2.2. The single Antigen case By using a list for each cellular entity, we prevent the exponential growth of memory requirements for cellular data structures but, unfortunately, the problem still remains for the molecules. It would be useless to assign a unique identity to each Ab to distinguish between two functionally identical molecules. Moreover, there are many more molecules (Ab, Ag and immuno complexes) than cells. Finally, the CS automaton, like other models [6], allows to model mutating Ag (like AIDS virus). This entails that the code must store the information required to simulate all possible run-time mutations of the Ags. So, the computational space complexity still remains exponential with the bit string length `, as far as the Abs and Ag are concerned. However, at least so far, we have been interested in observing the immune response just against predetermined Ags (which do not mutate). In this case it is possible to reduce dramatically the complexity of the code. Instead of considering the whole range of possible values for the Abs strings, a Hash table is exploited to pick out, from the data structures representing the stimulated PLB cells and the Ab, just those entries needed for the response to the Ag. This is useful because for matching values lower than mc there is no stimulation. Since the MinMatch is always greater than `/2, the number of possibilities for the stimulation and proliferation of PLBs and

338

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

  ` ` Abs is equal to 6k=m that is significantly smaller than 2` . For the ` = 14, mc = 10 case, the complexity is c k reduced to ∼ 9% of the most general case. By using this technique, the parallel version of the simulator (ParImm) can handle millions of cells and molecules regardless of the grid size.

3. The parallel version of the simulator The parallel version of ImmSim has been developed in two phases. In the first and preliminary version of the code (ParImm1.0) we adopted a simple master–worker scheme. This was good enough from the performance viewpoint but its storage requirements prevented us from running large simulations, namely the very motivation behind the development of a parallel version. Nevertheless it is useful to report here some details about ParImm1.0 to allow a better understanding of the choices we made for ParImm2.0. 3.1. Parallel Immsim1.0 In the serial ImmSim most of the time is spent in computing the interaction among entities belonging to the same site (insite). For ParImm1.0 we decided to run in parallel just that part, so there was a master and a set of workers. At each iteration: (i) the master distributed a subset of the sites M1 , . . . , MNP−1 to each worker W1 , . . . , WNP−1 ; (ii) master and workers carried out the simulation of the interaction among entities on their subset of sites; (iii) the workers sent back the results to the master; (iv) the master carried out all the other parts of the simulation (i.e. cells propagation, cells birth and death and so on). During the execution of Step 2 the tasks did not need to exchange any information. Actually the workers never talked each other. All the communication was based on two basic collective operations: (i) a scatter to distribute data (Step 1); (ii) a gather to collect them at the end of Step 3. Such communication pattern had the advantage of minimizing the changes with respect to the serial version of the code. Moreover, the tasks executed just two, although pretty big (about 60KB per site), message-passing operations per iteration. The communication phase was more bandwidth than latency bound and this was positive since our platform (IBM SP2) is very suitable for such codes. However, ParImm1.0 shows all its serious limitations as soon as the issues described in Section 2.2 are considered. There we showed how the size of the huge data structure which represent the molecules grows exponentially with respect to (`). In ParImm1.0 that static multidimensional array was not partitioned among the NP processors, as a consequence simulations with longer ` did not show any real advantage running in parallel since the memory swapping activity of the master was and remained the actual bottleneck. So there was a compelling need to split the data structures among all the processors in order to avoid the memory swapping on the master node and allow, in such way, to run more realistic (i.e. with a longer bit string length) simulations. 3.2. Parallel Immsim2.0 In the new version of the parallel code (ParImm2.0) all the phases of the simulation run in parallel. This major enhancement has been made possible by using a particular technique for the Domain Decomposition of the cellular automaton space. The basic idea is to consider each subdomain as an independent entity which interacts with other independent entities. In other words there is no single list split among the processors but as many lists as the number of processors in use.

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

339

The body is mapped onto a bi-dimensional grid of dimX × dimY sites, with periodic boundary conditions in both directions (up–down, left–right). However, to make the internal management of the lists easier, the sites of the automaton are not arranged as a two-dimensional array but as a linear chain. The transformation is carried out by a simple function (X, Y ) → N and does not change the global toroidal topology of the body. The assignment of the sites determines the ranking of the PEs with the exception of the first PE which always takes the first slice of the rectangular grid, (i.e. the sites from 0 to nsites − 1). In general the processor that works on the sites within the range [k · nsites, (k + 1) · · · nsites − 1] is indicated as PEk . Each PE manages private lists of the B, PLB, T and AP cells belonging to its subdomain and the data structures describing its molecules as if they were the whole body. No PE keeps a copy of lists or data structures belonging to other PEs and thus for a fixed ` the memory required on a PE decreases linearly as NP grows. Elements are deleted from or inserted in the lists of cells (and molecule arrays) when a subdomain boundary is crossed during the phase of diffusion. When a cell leaves the domain of a PE to migrate to one of the two nearest neighbour PEs, the corresponding element is deleted from the original list and all the relevant attributes are packed in a message sent to the PE which now owns the cell (PEd ). PEd unpacks the message and inserts the attributes of the incoming cell in a new element that becomes the head of the list. Actually each PE packs all its outgoing cells in just two messages for performance reasons. The first is directed to (PE − 1) mod NP, whereas the second goes to (PE + 1) mod NP. 3.3. ParImm2.0 results Both ParImm1.0 and ParImm2.0 have been realized with the primitives for parallel programming defined by the PVM software package [7]. The main advantage of PVM is, for our purpose, the very simple mechanism for packing multiple and heterogeneous elements in a single message. The results reported here have been produced on an IBM SP2 equipped with the high performance switch (HPS) and a proprietary version of PVM named PVMe [8]. All the timings we report are the total elapsed time of the simulation. For a fair comparison the single processor timings are obtained with the serial version of the simulator. A number of runs with different grid sizes and (`) have been executed. In Fig. 4 we present the results for two grids (16 × 16 and 16 × 16) and three values of ` (8, 10 and 12) on 1, 2, 4, 8 and 16 PEs. To keep the concentration of cells constant in a single site the total initial number of entities is doubled going from the small (16 × 16) to the large (16 × 32) test case.

Fig. 4. Time and SpeedUp (T (1)/T (NP)) for several test cases.

340

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

Fig. 5. B-cell and Antibodies population for a large test case (grid size 16 × 32, ` = 14).

From Fig. 4 it is apparent that the execution time increases linearly with respect to the size of the automaton whereas a longer ` determines an exponential growth. We knew that the average message size (140KB) was large enough to absorb the detrimental effects of the latency so we forecasted a good communication performance. Instrumented runs confirmed that the rate of the messages exceeds 30MB/s. Such communication rate along with a favourable ratio between computation and communication explains the fairly good speed-up observed, in particular for the largest simulations (` = 12). Since in ParImm2.0 all the parts of the simulation run in parallel, good scalability should be achievable on any number of PE provided that ` is long enough. A look at the results produced with a larger grid (16 × 32) and ` = 14 (Fig. 5) confirms the behaviour described in Section 2.1. However, the higher specificity (214 ) allows to gain a better insight into phenomena like the hypermutation of Abs [9]. 3.4. Porting to other platforms The choice of the PVM programing model has made the porting of ParImm to other platforms fairly easy. Indeed, besides the original SP2 version, ParImm runs on a number of other platforms: from network of workstations (NOWs) to the Cray/T3E supercomputer. However, the public domain version of PVM has been used just on clusters. For performance reasons, for instance, we have been forced to use our own version of PVM on the Sun Ultra Sparc. This version is similar to PVMe in design but it is specially tuned for shared memory architectures [10]. On the T3E we have resorted to the Cray PVM environment. The only problem with the Cray PVM is that it follows strictly the SPMD philosophy, the pvm spawn() primitive is not supported and there is a one-to-one correspondence between task and processor. These features have required slight changes to the ParImm code. We expect the T3E to be our main platform for ‘large’ simulations (here large means ` greater than 20 and number of cells in the order of millions).

4. Future work Cellular automata are usually understood as discrete dynamical systems governed by uniform evolution laws, i.e. all sites move in time according to the same update rules (lockstep mode). With ParImm, it is possible to have different interaction rules on each site. This is potentially useful in view of future applications in which a non-uniform automaton would be used to model distinct parts of the IS.

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

341

Our highly tuned implementation of the Celada–Seiden model should allow us to run simulations with millions of cells and a degree of specificity close to the real value. Besides the obvious interest from the biological viewpoint, the study of the statistical properties of such large dynamical systems bears an interest on its own. In a recent paper, we have analyzed the learning process of immune system in terms of statistical observables such as the entropy and the information gain [11] in a relatively small system ` = 12. The results look promising, in the sense that these observables seem to capture some degree of universality associated with collective modes of the IS response. Of course, much more numerical evidence needs to be produced in the future to back up such a bold statement. In particular, the robustness of the aforementioned collective modes needs to be tested over a broader range of parameters, such as the grid size, number of cells and specificity repertoire. We believe our parallel version of the Celada–Seiden model is ideally placed to achieve this goal.

Acknowledgements F. Celada and P. Seiden are kindly acknowledged for useful discussions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

D. Stauffer, R. Pandey, Comput. Phys. 6(4) (1992) 404. P. Seiden, F. Celada, J. Theor. Biol. 158 (1992) 329–357. F. Celada, P. Seiden, Immunol. Today 13(2) (1992) 56–62. M. Kaufman, J. Urbain, R. Thomas, J. Theor. Biol. 114 (1985) 527. R. Pandey, D. Stauffer, J. Statis. Phys. 61 (1990) 235. R.M. Zorzenon dos Santos, A.T. Bernades, Immunization and aging: a learning process in the immune network, Physical Rev. Lett. 81 (14) (1998). V.S. Sunderam, G.A. Geist, J. Dongarra, R. Manchek, Parallel Computing 20 (4) (1994). M. Bernaschi, G. Richelli, Proc. HPCN Europe 1995, Bob Hertzberge, Giuseppe Serazzi (Eds.), Lecture Notes in Computer Science, Springer, Berlin, p.919. F. Celada, P.E. Seiden, Eur. J. Immunol. 26 (1996) 1350–1358. M. Bernaschi, Efficient message passing on UNIX shared memory multiprocessors, Future Gen. Comput. Sys. J., accepted for publication. M. Bernaschi, F. Castiglione, S. Succi, Collective dynamics in the immune system response, Phys. Rev. Lett. 79 (22) (1997).

Massimo Bernaschi graduated in physics in 1987 at ‘Tor Vergata’ University in Rome. During and after the thesis work he took part in a project of the National Institute for Nuclear Physics (called APE collaboration) aimed at developing a parallel computer for lattice gauge theory. At the end of 1987 he joined the IBM European Center for Scientific and Engineering Computing (ECSEC) in Rome. He has spent ten years with IBM working with a number of hardware and software platforms for parallel and distributed computing. Currently he is with the Italian National Research Council (CNR) as a director of technology office in the Institute for Computing Applications.

Filippo Castiglione graduated in computer science in 1993 at the University of Milan. Besides his work on the simulation of the immune system, he works on projects involving Image Analysis and Pattern Recognition in biomedical science. He is currently a PhD student at the Center for Parallel Computing at the University of Cologne, Germany.

342

M. Bernaschi et al. / Future Generation Computer Systems 15 (1999) 333–342

Sauro Succi graduated in nuclear engineering at the University of Bologna and holds a PhD in plasma physics from the Swiss Federal Institute of Technology in 1987. In 1986 he joined IBM ECSEC (European Center for Scientific and Engineering Computing), where he served as Research Scientist, Senior Research Scientist and Industry Sector Coordinator. In 1995, he has been appointed Research Director at the Institute for Computing Applications of the Italian National Research Council, where he is currently leading a group in applied mathematics. Dr. Succi has published extensively in plasma physics, kinetic theory, computational fluid dynamics and parallel computing. He has been holding visiting positions at the University of Paris, Chicago and MIT Cambridge. He serves as an associate editor of the International Journal of Modern Physics C: Physics and Computers and Journal of Scientific Computing. He’s member of the American Society for Industrial and Applied Mathematics, the European Physical Society, the American Physical Society and the New York Academy of Sciences.