A Multiscale Agent-based Model of Morphogenesis in Biological Systems

Abstract

-Studying the complex phenomenon of pattern formation created by the gene expression is a big challenge in the field of developmental biology. This spatial self-organisation autonomously emerges from the morphogenetic processes and the hierarchical organisation of biological systems seems to play a crucial role. Being able to reproduce the systems dynamics at different levels of such a hierarchy might be very useful. In this paper we propose the adoption of the agent-based model as an approach capable of capture multi-level dynamics. Each cell is modelled as an agent that absorbs and releases substances, divides, moves and autonomously regulates its gene expression. As a case study we present an agent-based model of Drosophila melanogaster morphogenesis. We then propose a formalisation of the model which clearly describe its main components. We finally show simulation results demonstrating the capability of the model of reproducing the expression pattern of the embryo.

I. INTRODUCTION

Developmental biology is an interesting branch of life science that studies the process by which organisms develop, focussing on the genetic control of cell growth, differentiation and movement. A main problem in developmental biology is understanding the mechanisms that make the process of vertebrates’ embryo regionalisation so robust, making it possible that from one cell (the zygote) the organism evolves acquiring the same morphologies each time. This phenomenon involves at the same time the dynamics of – at least – two levels, including both cell-to-cell communication and intracellular phenomena: they work together, and influence each other in the formation of complex and elaborate patterns that are peculiar to the individual phenotype. This happens according to the principles of downward and upward causation, where the behaviour of the parts (down) is determined by the behaviour of the whole (up), and the emergent behaviour of the whole is determined by the behaviour of the parts [ 19 ].

Modelling embryo- and morphogenesis presents big challenges: (i) there is a lack of biological understanding of how intracellular networks affect multicellular development and of rigourous methods for simplifying the correspondent biological complexity: this makes the definition of the model a very hard task; (ii) there is a significant lack of multilevel models of vertebrate development that capture spatial and temporal cell differentiation and the consequent heterogeneity in these four dimensions; (iii) on the computational framework side, there is the need of tools able to integrate and simulate dynamics at different hierarchical levels and spatial and temporal scales.

A central challenge in the field of developmental biology is to understand how mechanisms at intracellular and cellular level of the biological hierarchy interact to produce higher level phenomena, such as precise and robust patterns of gene expressions which clearly appear in the first stages of morphogenesis and develop later into different organs. How does local interaction among cells and inside cells give rise to the emergent self-organised patterns that are observable at the system level?

The above issues have already been addressed with different approaches, including mathematical and computational ones. Mathematical models, on the one side, are continuous, and use differential equations—in particular, partial differential equations describing how the concentration of molecules varies in time and space. A main example is the reaction-diffusion model developed by [ 18 ] and applied to the Drosophila melanogaster (Drosophila in short) development by [ 15 ]. The main drawback of mathematical models is the inability of building multi-level models that could reproduce dynamics at different levels.

Computational models, on the other side, are discrete, and model individual entities of the system—cells, proteins, genes. The agent-based approach is an example of such a kind of models. Agent-based modelling (ABM) is a computational approach that can be used to explicitly model a set of entities with a complex internal behaviour and which interact with the others and with the environment generating an emergent behaviour representing the system dynamics. Some work has already been done which applies ABM in morphogenesislike scenarios: a good review is proposed in [ 17 ]. Most of these models generate artificial pattern – French and Japanese flags [ 2 ] – realising bio-inspired models of multicellular development in order to obtain predefined spatial structures. At the best of our knowledge, however, few results have been obtained till now in the application of ABM for analysing real phenomena of morphogenesis.

In order to get the benefits of both approaches, hybrid frameworks has been developed. For instance, COMPUCELL 3D [ 3 ] combines discrete methods based on cellularautomata to model cell interactions and continuous model based on reaction-diffusion equation to model chemical diffusion. COMPUCELL 3D looks like a very promising framework whose main limitation is represented by the lack of a suitable model for cell internal behaviour—gene regulatory network in particular.

In this paper we present an agent-based model of the Drosophila embryo development, reproducing the gene regulatory network that causes the early (stripes-like) regionalisation of gene expression in the anteroposterior axis [ 21 ], [ 15 ]. The embryo is modelled as a set of agents, where each agent is a cell. Our approach allows the gene-regulatory network to be directly modelled as the internal behaviour of an agent, whose state reproduces the gene expression level and dynamically changes according to functions that implement the interactions among genes. It also allows the cell interacting capability mediated by morphogens to be modelled as the exchange of messages among agents that absorb and secrete – from and towards the environment – the molecules that are then able to diffuse over the environment.

The remainder of this paper is organised as follows: The role of hierarchy in the spatial self-organisation of gene expression during morphogenesis is first highlighted along with the main biochemical mechanisms taking place in this phenomenon. The agent-based approach is then presented with the modelling abstractions it provides. The third part describes the biological principles of Drosophila embryo development, while the fourth part reports the ABM we have developed, formalised and implemented. Simulation results are then discussed, followed by concluding remarks.

II. THE ROLE OF HIERARCHY IN MORPHOGENESIS Complex systems in general exhibit a hierarchical organisation that divide the system into levels composed by many interacting elements whose behaviour is not rigid, and is instead self-organised according to a continuous feedback between levels. Hierarchy has therefore a crucial role in the static and dynamic characteristics of the systems themselves. An example is given by biological systems: an outstanding property of all life is the tendency to form multi-levelled structures of systems within systems. Each of these forms a whole with respect to its parts, while at the same time being a part of a larger whole. Biological systems have different level of hierarchical organisation – (1) sequences; (2) molecules; (3) pathways (such as metabolic or signalling); (4) networks, collections of cross-interacting pathways; (5) cells; (6) tissues; (7) organs – and the constant interplay among these levels gives rise to their observed behaviour and structure. This interplay extends from the events that happen very slowly on a global scale right down to the most rapid events observed on a microscopic scale. A unique molecular event, like a mutation occurring in particularly fortuitous circumstances, can be amplified to the extent that it changes the course of evolution. In addition, all processes at the lower level of this hierarchy are restrained by and act in conformity to the laws of the higher level.

In this contest, an emblematic process is morphogenesis, which takes place at the beginning of the animal life and is responsible for the formation of the animal structure. Morphogenesis phenomena includes both cell-to-cell communication and intracellular dynamics: they work together, and influence each other in the formation of complex and elaborate patterns that are peculiar to the individual phenotype.

A. The biology of development

Animal life begins with the fertilisation of one egg. During the development, this cell undergoes mitotic division and cellular differentiation to produce many different cells. Each cell of an organism normally owns an identical genome; the differentiation among cells is then not due to different genetic information, but to a diverse gene expression in each cell. The set of genes expressed in a cell controls cell proliferation, specialisation, interactions and movement, and it hence corresponds to a specific cell behaviour and role in the entire embryo development.

One possible way for creating cells diversity during embryogenesis is to expose them to different environmental conditions, normally generated by signals from other cells, either by cell-to-cell contact, or mediated by cues that travel in the environment.

On the side of intracellular dynamics, signalling pathways and gene regulatory networks are the means to achieve cells diversity. Signalling pathways are the ways through which an external signal is converted into an information travelling inside the cell and, in most of the cases, affecting the expression of one or more target genes. The signalling pathways are activated as a consequence of the binding between (i) a cue in the environment and a receptor in the cell membrane, or (ii) two membrane proteins belonging to different cells. The binding causes the activation of the downstream proteins until a transcription factor that activates or inhibits the expression of target genes.

During morphogenesis few pathways are active. They work either as mutual inhibitors, or as mutual enhancers. The idea is that there are regions where the mutual enhancers are active and interact giving rise to positive feedbacks. Pathways active in different regions work probably as mutual inhibitors. There are then boundary regions where we can observe a gradient of activity of the different sets of pathways, due to the inhibitory effect of the pathways belonging to neighbour regions.

III. THE AGENT-BASED APPROACH

In literature, agent-based systems – in particular MultiAgent Systems (MAS) – are considered as an effective paradigm for modelling, understanding, and engineering complex systems, providing a basic set of high level abstractions that makes it possible to directly capture and represent the main aspects of such complex systems, such as interaction, multiplicity and decentralisation of control, openness and dynamism [ 13 ], [ 12 ], [ 11 ]. A MAS can be characterised by three key abstractions: agents, societies and environment. Agents are the basic active components of the systems, executing pro-actively and autonomously. Societies are formed by set of agents that interact and communicate with each other, exploiting and affecting the environment where they are situated. Such an environment plays a fundamental role, as a context enabling, mediating and constraining agent activities [ 20 ].

By adopting an agent-based approach, biological systems can be modelled as a set of interacting autonomous components – i.e., as a set of agents –, whereas their chemical environment can be modelled by suitable agent environment abstractions, enabling and mediating agent interactions. In particular, MAS provide a direct way to model: (i) the individual structures and behaviours of different entities of the biological system as different agents (heterogeneity); (ii) the heterogeneous – in space and time – environment structure and its dynamics; (ii) the local interactions between biological entities/agents (locality) and their environment. An agentbased simulation means executing the MAS and studying its evolution through time, in particular: (i) observing individual and environment evolution; (ii) observing global system properties as emergent properties from agent-environment and inter-agent local interaction; (iii) performing in-silico experiments. The approach is ideal then for studying the systemic and emergent properties that characterise a biological system, which are meant to be reproduced in virtuo. In the context of biological system, agent-based models can therefore account for individual cell biochemical mechanisms – gene regulatory network, protein synthesis, secretion and absorption, mitosis and so on – as well as the extracellular matrix dynamic – diffusion of morphogens, degradation and so on – and their dynamic influences on cell behaviour.

IV. THE DROSOPHILA EMBRYO DEVELOPMENT One of the best example of pattern formation during morphogenesis is given by the patterning along the anteroposterior axis of the fruit fly Drosophila melanogaster.

A. Biological background

The egg of Drosophila is about 0.5 mm long and 0.15 mm in diameter. It is already polarised by differently localised mRNA molecules which are called maternal effects The early nuclear divisions are synchronous and fast (about every 8 minutes): the first nine divisions generate a set of nuclei, most of which move from the middle of the egg towards the surface, where they form a monolayer called syncytial blastoderm. After other four nuclear divisions, plasma membranes grow to enclose each nucleus, converting the syncytial blastoderm into a cellular blastoderm consisting of about 6000 separate cells.

Up to the cellular blastoderm stage, development depends largely – although not exclusively – on maternal mRNAs and proteins that are deposited in the egg before fertilisation. After cellularisation, cell division continues asynchronously and at a slower rate, and the transcription increases dramatically. Once cellularisation is completed the gene expression regionalisation is already observable.

The building blocks of anterior-posterior axis patterning are laid out during egg formation thanks to the maternal effects. Bicoid and caudal are the maternal effect genes that are most important for patterning of anterior parts of the embryo in this early stage. They are transcription factors that drive the expression of gap genes such as hunchback (Hb), Kru¨ppel (Kr), knirps (Kni) and giant (Gt), as shown in the diagram of Fig. 1; there, tailless (Tll) also appears as gap genes whose regulation we do not represent here. Gap genes together with maternal factors then regulate the expression of downstream targets, such as the pair-rule and segment polarity genes. The segmentation genes specify 14 parasegments that are closely related to the final anatomical segments [ 1 ], [ 6 ].

Our model consists of a set of agents that represent the cells, as well as of a grid-like environment representing the extracellular matrix. Agent internal behaviour reproduces the gene regulatory network of the cell, while agent interaction with the environment models the process of cell-to-cell communication mediated by the signalling molecules secreted in and absorbed by the extra-cellular matrix. Our model aims at reproducing the expression pattern of the gap genes, before the pair-rule genes are activated.

A. Model of the cell

We model different cell processes: secretion-absorption diffusion of chemicals from and towards the environment, cell growth, cell movement and cell internal dynamics—gene regulatory network in particular.

1) Chemical diffusion: We model the process of molecule secretion and absorption as facilitated diffusion—the literature lacks of information about the transport mechanisms of such transcription factors and about the rate of diffusion.

2) Gene regulatory network: Gene transcription begins with the binding at the gene promoter of one or more transcription factors. Gene transcription might also be repressed once transcription factors bind to other control regions called silencers. This activation/inhibition is stochastic [ 10 ] and highly depends on the concentration of transcription factors. For those genes whose transcription is regulated by a set of other gene products we define a probability of transcription as a sum of positive and negative contributions from the concentration of enhancers and silencers, respectively. For instance, the probability of transcription of hunckback, according to the graph of Fig. 1, is then calculated as:

Ph = f ([Bicoid ]) + f ([Hunchback ]) + f ([Tailless ]) f ([Knirps ]) f ([Kruppel ]) where f is a linear function with the proportionality constant representing the strength of interaction. Then if Ph > 0 the protein is synthesised, otherwise the gene remains silent.

No distinction has been done in the model between anterior (a) and posterior (p) hunckback and giant, whose different expression only deals with the spatial distribution of maternal products.

3) Movement: The model of cell movement considers the chemotaxis phenomenon which is known to be responsible for cell sorting during morphogenesis [ 4 ]. This model component is inspired at a previous work that considers chemotaxis as an important actor for the creation of self-organised structures [ 5 ]. We here assume that cells move towards the direction created by the gradient of the protein which is more expressed in each of them.

4) Mitosis: According to Fig. 2 where we show how the number of cells varies in the first four hours of embryo development – until the cleavage cycle 14, temporal class 8 – we computed the rate of division as a function of time: cell division is fast and synchronous until cleavage cycle 9, then slows down and becomes asynchronous. The rate of division is constant in the first hours of development (9.05 min 1), then decreases until a low value (0.2 min 1), as it appears in Figure 3.

6000 4000 2000 1000 lls500 e c f ro100 be 50 m u N 1 2 3 4 5 6 7 8 9 10 11 12 13 14A−114A−214A−134A−414A−5 14A−6

14A−8 14A−7

B. Model of the environment

The 3D-tapered structure of the embryo, as in Figure 4, is modelled as a 2D-section of the embryo along the anteroposterior axis (c) under the assumption that the dynamics along the other two axis, a and b, does not influence what happens along the c axis. The space scale is 1:3.33 according to the real dimension of the embryo where the antero-posterior axis is almost three times the dorso-ventral one a. Space is not continuous but grid like, and each location might be occupied both by a set of morphogenes and by a cell.

The environment has its own dynamics, which mainly consists in the diffusion of morphogenes from region with bigger 20 40 60 80 concentration to region with lower concentration, according to the Fick’s low that the diffusive flux is proportional to the local concentration gradient [ 16 ]. This law is used in its discretised form.

C. Model formalisation

Figure 5 shows some of the statecharts [ 9 ] used to formally describe the cell behaviour, in term of a hierarchical structure of states and event-triggered transitions. The main macro state Alive contains two states, Placing and Life Cycle. As soon as a cell is created (Alive macro state), first it moves to find its place inside the embryo (Placing sub-state) and then – when the moving is completed – it starts a life-cycle (Life Cycle sub-state) until its death. Such a life cycle is modelled with two parallel processes that concern proteins’ expulsion and absorbing – composed by Rest A (idle), Expelling and Absorbing sub-states – and cell protein synthesis and cell mitosis – composed by sub-states Rest B (idle), Synthesizing and Mitosis Cycle sub-states. The cell expels or absorbs molecules from its environment depending on the concentration of such substances inside the cell and the one perceived in its local environment. The protein synthesis process is triggered by the perception of a specific activating protein (evActiveProtein event in the statechart), while mitosis is triggered by a specific situation which includes timing and other cell state conditions (evMitosis).

D. Model implementation and simulation procedure

The model is implemented on top of Repast Simphony1, an open-source, agent-based modelling and simulation toolkit. 1http://repast.sourceforge.net/index.html It provides all the abstraction for directly modelling the agent behaviour and the environment. It implements a multithreaded discrete event scheduler.

Simulations are executed from the cleavage cycle 11, when the zygotic expression begins. We used the experimental data available online in the FlyEx database2. The data contains quantitative wild-type concentration profiles for the protein products of the seven genes – Bcd, Cad, Hb, Kr, Kni, Gt, Tll – during cleavage cycles 11 up to 14A, which constitutes the blastoderm stage of Drosophila development. These data are used to validate the model dynamic. Expression data from cleavage cycle 11 are used as initial condition—see Fig. 7. The concentration of proteins are unitless, ranging from 0 to 255, at space point x, ranging from 0 to 100 % of embryo length.

Model parameters are: (i) diffusion constants of morphogenes motion; (ii) rates of gene interactions; (iii) rates of protein synthesis. Few data are available in literature for inferring the diffusion constants. We took inspiration from the work of [ 7 ] that calculates the diffusion rate for Bicoid and we imposed the value for all the morphogenes at 0.3 m2=sec. The rates of gene interactions and of protein synthesis are determined through a process of automatic parameter tuning. The task is defined as an optimisation problem over the parameter space. The optimisation makes use of metaheuristics – particle swarm optimisation – to find a parameter configuration such that the simulated system has a behaviour comparable with the real one [ 14 ]. We supported the automatic parameter tuning with a process of model refinement which slightly changed 2http://flyex.ams.sunysb.edu/flyex/index.jsp the topology of gene regulatory network, adding some edges that we found necessary for obtaining the real behaviour. An argumentation about the final model is provided in the Discussion.

Qualitative results charted in the 2D-grid are shown in Fig. 6 (top) for expression of hb, kni, gt, Kr at the eighth time step of cleavage cycle 14A. The image shows for each cell of the embryo the genes with higher expression. It clearly displays the formation of a precise spatial pattern along the A-P axis but it does not give any information about gene expression level. Experimental data are also provided in Fig. 6 (bottom) with 2D-Atlas reconstructing the expression level of the four genes in A-P sections of the embryo. More precise information about simulation behaviour are given with the quantitative results 10 20 30 40 60 70 80 90 100

00 l200 e v len150 o i s rse100 p x e ne 50 e G 00 l200 e v len150 o i s rse100 p x e ne 50 e G 00 10 20 30 40 50 60 70 80 90 100 provided in Fig. 8. A comparison shows that the expression pattern of genes Hb, Kni, Gt and Kr nicely fit the spatial distribution shown in the experimental data: Hb is expressed in the left pole until about 45% of embryo length, while it does not appear on the right as it should between about 85% and 95%; Kni is correctly expressed on the extreme left and between 65% and 75% but it is slightly over-expressed on the right; Gt is reproduced in the correct regions but overexpressed in the extreme left and slightly under-expressed between 20% and 30%; finally, Kr properly appears between 40% and 60%.

VI. DISCUSSION

Through the model refinement we found the network showed in Fig. 9 where some more interactions are performed. The weight in sec 1 of each node is then reported in Fig. 10.

Bcd and Cad are activators of the gap genes. As maternal factor their central role is in fact to input the wave of zygotic expression. In particular, given the spatial distribution of their expression, Bcd is responsible for the activations on the left side of the embryo, while Cad in the opposite side. Tll enhances Hb expression while inhibits the expression of all the others as in the previous model. The interactions among gap genes are slightly different. As before Hb and Kni on one side and Gt and Kr on the other side inhibits one each other, and from the parameters found we infer that these are the strongest inhibitions among gap genes; Hb then weakly inhibits Kr and vice-versa, as well as Gt versus Kni. New weak edges have been found between Kni versus Gt, and Kr versus Kni.

As far as we know, there are no evidences in biological literature that already support the above results. It might be a starting point for new laboratory experiments.

VII. CONCLUSION

The process of spatial organisation resulting from the morphogenesis process is demonstrated to be highly-dependent by the interplay between the dynamics at different levels of the biological systems hierarchical organisation. In modelling and simulating the phenomena of morphogenesis it might be appropriate to reproduce such a hierarchy. In this work we have described the application of ABM as an approach capable of supporting multi-level dynamics.

We studied the phenomenon of pattern formation during Drosophila embryo development, modelling the interactions between maternal factors and gap genes that originate the early regionalisation of the embryo. The possibility to model both the reactions taking place inside the cells that regulate the gene expressions, and the molecules diffusion that mediates the cellto-cell communication, makes it possible the reproduction of the interplay between the two levels in order to verify its fundamental role in the spatial self-organisation characteristic of such a phenomenon.

The model is formally described using the statecharts that make it possible to clearly show the model components and how they behave.

The simulation results presented show the formation of a precise spatial pattern which have been successfully compared with observations acquired from the real embryo gene expressions.

Future work will be firstly devoted to extending the model with the introduction of new phenomena on the side of both intracellular dynamics and cell-to-cell interaction. Gene regulatory network will be enlarged with other sets of genes which are downstream to gap genes such as the pair rule genes, even-skipped as first, whose expression gives rise at the characteristic segments of Drosophila embryo. Other mechanism regulating cell movements will then be added – for instance cell adhesion and mechanic forces – as soon as they are known to play a crucial role in cell sorting during morphogenesis.

Finally. we are planning to exploit the predictive power of the model analysing embryos that are not wild type, for instance performing in-silico Knock-Out experiments.

References

[1] B. Alberts , A. Johnson , J. Lewis , M. Raff , K. Roberts , and P. Walter . Molecular Biology of the Cell. Garland Science Textbooks. Garland Science, 4th edition , June 2002 .

[2] G. Beurier , F. Michel , and J. Ferber . A morphogenesis model for multiagent embryogeny . In L. M. Rocha , L. S. Yaeger , M. A. Bedau , D. Floreano , R. L. Goldstone , and A . Vespignani, editors, Artificial Life X: Proceedings of the Tenth International Conference on the Simulation and Synthesis of Living Systems , pages 84 - 90 . MIT Press, Cambridge, MA, 2006 .

[3] T. M. Cickovski , C. Huang , R. Chaturvedi , T. Glimm , H. G. E. Hentschel , M. S. Alber , J. A. Glazier , S. A. Newman , and J. A. Izaguirre . A framework for three-dimensional simulation of morphogenesis . IEEE/ACM Transactions on Computational Biology and Bioinformatics , 2 : 273 - 288 , 2005 .

[4] J. Davies . Mechanisms of Morphogenesis. Academic Press, Oct. 2005 .

[5] M. Eyiyurekli , L. Bai , P. I. Lelkes , and D. E. Breen . Chemotaxis-based sorting of self-organizing heterotypic agents . In SAC '10: Proceedings of the 2010 ACM Symposium on Applied Computing , pages 1315 - 1322 , New York, NY, USA, 2010 . ACM.

[6] S. F. Gilbert . Developmental Biology, Eighth Edition . Sinauer Associates Inc., Eighth edition , Mar. 2006 .

[7] T. Gregor , E. Wieschaus , A. McGregor , W. Bialek , and D. Tank . Stability and nuclear dynamics of the bicoid morphogen gradient . Cell , 130 ( 1 ): 141 - 152 , July 2007 .

[8] V. V. Gursky , J. Jaeger , K. N. Kozlov , J. Reinitz , and A. M. Samsonov . Pattern formation and nuclear divisions are uncoupled in drosophila segmentation: comparison of spatially discrete and continuous models . Physica D: Nonlinear Phenomena , 197 ( 3-4 ): 286 - 302 , October 2004 .

[9] D. Harel . Statecharts: A visual formalism for complex systems . Sci. Comput. Program. , 8 ( 3 ): 231 - 274 , 1987 .

[10] M. Kaern , T. C. Elston , W. J. Blake , and J. J. Collins. Stochasticity in gene expression: from theories to phenotypes . Nature reviews. Genetics , 6 ( 6 ): 451 - 464 , June 2005 .

[11] F. Klu ¨gl, C. Oechslein , F. Puppe , and A. Dornhaus . Multi-agent modelling in comparison to standard modelling . In F. J. Barros and N. Giambiasi, editors, Artificial Intelligence, Simulation and Planning in High Autonomy Systems , pages 105 - 110 . SCS Publishing House, 2002 .

[12] E. Merelli , G. Armano, N. Cannata , F. Corradini , M. d'Inverno , A. Doms , P. W. Lord , A. Martin , L. Milanesi , S. Mo¨ller, M. Schroeder, and M. Luck . Agents in bioinformatics, computational and systems biology . Briefings in Bioinformatics, 8 ( 1 ): 45 - 59 , 2007 .

[13] F. Michel , J. Ferber , and A. Drogoul . Multi-Agent Systems and Simulation: a Survey From the Agents Community's Perspective . In A. U. Danny Weyns, editor, Multi-Agent Systems: Simulation and Applications , Computational Analysis, Synthesis, and Design of Dynamic Systems , pages 3 - 51 . CRC Press - Taylor & Francis, 05 2009 .

[14] S. Montagna and A. Roli . Parameter tuning of a stochastic biological simulator by metaheuristics . In AI*IA , volume 5883 of Lecture Notes in Computer Science, pages 466 - 475 . Springer, 2009 .

[15] T. J. Perkins , J. Jaeger , J. Reinitz , and L. Glass . Reverse engineering the gap gene network of Drosophila Melanogaster . PLoS Comput Biol , 2 ( 5 ):e51, 05 2006 .

[16] W. Smith and J. Hashemi . Foundations of Materials Science and Engineering . McGraw-Hill , 4th edition , July 2005 .

[17] B. Thorne , A. Bailey , D. Desimone , and S. Peirce . Agent-based modeling of multicell morphogenic processes during development . Birth Defects Res C Embryo Today , 81 ( 4 ): 344 - 353 , 2008 .

[18] A. M. Turing . The chemical basis of morphogenesis . Philosophical Transactions of the Royal Society (B) , 237 : 37 - 72 , 1952 .

[19] A. M. Uhrmacher , D. Degenring , and B. Zeigler . Discrete event multilevel models for systems biology . In C. Priami, editor, Transactions on Computational Systems Biology I, volume 3380 of Lecture Notes in Computer Science , pages 66 - 89 . Springer, 2005 .

[20] D. Weyns , A. Omicini , and J. Odell . Environment as a first-class abstraction in multi-agent systems . Autonomous Agents and Multi-Agent Systems , 14 ( 1 ): 5 - 30 , Feb. 2007 . Special Issue on Environments for Multi-agent Systems .

[21] D. Yamins and R. Nagpal . Automated global-to-local programming in 1-d spatial multi-agent systems . In 7th International Joint Conference on Agents and Multi-Agent Systems (AAMAS-08) , pages 615 - 622 , Estoril, Portugal, 12 -16May 2008 . IFAAMAS.