Agent-based modeling: case study in cleavage furrow models.

The number of studies in cell biology in which quantitative models accompany experiments has been growing steadily. Roughly, mathematical and computational techniques of these models can be classified as "differential equation based" (DE) or "agent based" (AB). Recently AB models have started to outnumber DE models, but understanding of AB philosophy and methodology is much less widespread than familiarity with DE techniques. Here we use the history of modeling a fundamental biological problem-positioning of the cleavage furrow in dividing cells-to explain how and why DE and AB models are used. We discuss differences, advantages, and shortcomings of these two approaches.

MATHEMATICAL AND COMPUTATIONAL MODELING IN CELL BIOLOGY

One of the main goals of cell biology is to understand how a cell works as a machine: where and how key proteins interact and achieve a desired function. It is no surprise, then, that an experimental article very often ends with a final figure that distills insights from results into a cartoon-like qualitative model—a mechanistic blueprint depicting protein machinery at work. More often than not, this is all that is needed, but there are cases in which there is a clear question that requires a mathematical or computational model. There could be various incentives for modeling—for example, to make sure that the cartoon does not contradict the rules of physics and chemistry, or the cartoon could be so complex that mathematics rather than limited intuition is necessary to demonstrate that the protein machine works as hypothesized (Mogilner ).

As experimental biology becomes more quantitative and complex, the number of cases in which modeling accompanies experimental studies grows. There are, however, many difficulties with using mathematics in biology. One of them is that modeling approaches have diversified so much recently that it is not easy to grasp how models work. If one looks at the greatest early successes of mathematical cell biology, it is hard to overlook two iconic examples: in the middle of the 20th century, Turing (1952) demonstrated mathematically that two chemicals, a slowly diffusing activator and a rapidly diffusing inhibitor, can generate spatially periodic patterns, which created a fundamental paradigm for morphogenesis. At around the same time, Hodgkin and Huxley (1952) showed numerically that a complex combination of ion fluxes through highly nonlinear voltage-gated channels can account for excitable electric pulses in nerve cell membrane.

These studies and hundreds of others used the powerful technique of differential equations (DEs), which goes back to Isaac Newton, who discovered that laws of mechanics can be expressed by relating rates of change of position and velocity (time derivatives) to position-dependent forces. Solutions of resulting ordinary differential equations (ODEs) turned out to be spectacularly effective in astronomy and physics. Later, Fourier proposed a diffusion equation—one of the most important partial differential equations (PDEs)—according to which the rate of change of a concentration at a given point is proportional to the spatial derivative of the spatial gradient of this concentration. DEs proved to be highly successful in the natural sciences and engineering. Indeed, Newton said that the laws of nature are written in the form of DEs, but he mostly meant physics.

DIFFERENTIAL EQUATIONS AND AGENT-BASED MODELS

Until recently, a majority of models in biology used DEs. The reason is very simple: a great number of biological problems deal with many molecules that move randomly and/or directionally and engage in chemical reactions (Figure 1). Very often, one can reasonably well describe the molecular ensemble by a continuous concentration/density function; the local rate of change of density can be described by a sum of a diffusion term (second spatial derivative of density, accounting for the random movements), a drift term (first spatial derivative of density, accounting for the directional movements), and an algebraic reaction term accounting for chemical reactions (Murray, 2002). A powerful analytical and numerical apparatus was developed over centuries to solve and understand DEs, and the comparison of the model’s predictions with experiments is in principle straightforward: DE solutions for the density as functions of time and the spatial coordinates could be compared with microscopy data.

FIGURE 1:

Comparison of DE and AB modeling. See the text for explanations.

Recently DE models started to lose out to so-called agent-based (AB) models. These are computational models simulating actions and interactions of autonomous agents (either individual molecules or collective molecular ensembles). The key notion on which AB models are based is that multiple agents interact according to simple rules; these rules are easy to encode in a computer program, and interactions are easy to simulate. However, despite their simplicity, these rules and interactions can generate complex emergent behavior.

To give an example of an AB model, one can consider the reaction-drift-diffusion system (see Figure 1 and the section Comparison of DE and AB Modeling). Whereas a DE model would deal with solving PDEs describing molecular density, an AB model would start with describing rules of behavior of individual molecules: at each small time interval, a molecule jumps randomly in space due to diffusion, shifts directionally due to drift, and, with a certain probability depending on the proximity of other molecules, disappears or appears or turns into another molecule. Then a great number of molecules would be simulated in parallel, giving us very complex spatial-temporal patterns (Figure 1 and Comparison of DE and AB Modeling).

The growing popularity of AB models is explained by many factors, one of which is that biological systems are, in fact, often characterized by complex emergent behavior, almost impossible to intuit, and based on simple interactions of a great number of molecular agents. Another reason is computers: relatively few DEs have analytical solutions. Turing, Hodgkin, and Huxley foresaw the paramount importance of computers for simulating complex sets of DEs that arise in biology. Advances in computing technology since then have made numerical solutions of DEs relatively easy and more recently have made possible challenging simulations of AB systems, finally bringing modeling to the experimental lab. Here we review the history of modeling of cleavage furrow positioning in cytokinesis to illustrate and compare DE and AB models as clearly, simply, and concisely as possible at the expense of not discussing subtleties and not being comprehensive. Other reviews fill this gap (Bryson ; An ; Holcombe ).

MODELING CLEAVAGE FURROW POSITIONING

Cell division at the end of mitosis starts with a cleavage furrow developing at the cell equator. Physically, the furrow is due to the contraction of the actomyosin ring at the cell cortex. Many serendipitous observations and ingenious micromanipulation experiments made over >100 years (Rappaport, 1961, 1971) have converged on four qualitative hypotheses-cartoons of the furrow positioning (Mishima, 2016; Figure 2A). Three of these models posit that the mitotic spindle activates, in a spatially precise way, the onset of the actomyosin ring assembly and contraction. The “astral stimulation model” (Rappaport, 1961) posits that astral microtubules (MTs) deliver a molecular activator to the cell actomyosin cortex; due to the overlap of the MT asters from the two sister centrosomes, the number of astral MTs at the cortex might peak at the equator, causing the furrow to form there. The “polar relaxation model” (Wolpert, 1960; White and Borisy, 1983), in contrast, proposes that the spindle restricts furrowing to the cell equator by means of the MT asters inhibiting cortical contractility at the cell poles. The “central spindle model” suggests that a molecular activator is delivered to the equator from a special bundle of MTs developing in the spindle midzone in anaphase. Finally, the “mitotic apparatus–independent model” hypothesizes that a mechanochemical pattern self-organizes in the cortex without the spindle as a central governor and positions the furrow to certain molecular cues in this pattern.

FIGURE 2:

(A) Four qualitative models of the furrow positioning. (B) Top, detailed astral stimulation model (Devore ); bottom, predicted distribution of the signaling molecule in the cell (modified from Devore ). (C) Snapshot from AB simulations in Odell and Foe (2008). MTs are green, centrosomes are red, and white dots show centralspindlin. Density plots surrounding the spherical cell show various averaged densities of centralspindlin (taken from Figure 1B in Odell and Foe, 2008).

A very clear quantitative question asking for mathematical modeling is immediately apparent when one looks at the cartoons of the astral stimulation and polar relaxation models: given certain geometries of the spindle MT asters and cell shape, would more signal be delivered to the equatorial or polar regions of the cortex? The former or latter could then be used as an argument in favor of the astral stimulation or polar relaxation model, respectively. This question prodded a few insightful DE modeling articles in the 1980s. White and Borisy (1983) computed the surface density of the astral MT plus ends at the cortex from two spindle asters in a spherical cell. In their simulations, the spindle was symmetrically positioned in the cell, and the MTs were stable and straight, stretching from the centrosome to the cortex. The strength of the signal from an individual MT was taken as proportional to an inverse power of the MT length (signal weakens with either centrosome–cortex distance, with MT splaying apart, or both). This study found that the maximal signal is at the poles, and therefore the polar relaxation mechanism is more likely.

However, the polar relaxation model has been shown to predict incorrectly the effects on cytokinesis of at least 15 experimental modifications of cellular shape (Devore ). For example, if the cell is an elongated cylinder with a small spindle in the middle, then clearly very little signal would be delivered to the cell poles. According to the polar relaxation model, this would imply furrows forming at the poles, in contradiction with experimental results. Devore therefore proposed the following modification based on visual examinations of spindle micrographs: there is a conical region for each aster in the middle of the spindle where the MTs do not grow, perhaps because they bump into the large structures of the opposite aster and destabilize (Figure 2B). In addition to this geometric assumption, Devore assumed that signaling molecules drift with a constant flux to the MT plus ends, detach from the plus ends near the cortex, diffuse, and are spontaneously degraded in the cytoplasm (Figure 2B). Then they translated these assumptions into a reaction-diffusion-drift PDE and solved this PDE numerically. The solutions of the PDE in such geometry predicted that the signal maximum at the poles of the spherical cell is only local; the global signal maxima are achieved at the equator (Figure 2B), leading the authors to conclude that the astral stimulation model accounts for the experimental data, whereas the data falsify the polar relaxation model. This and another modeling study (Harris and Harris, 1989) scanned a number of cell geometries and key model parameters (i.e., the angle θ defining the MT-void sector in the spindle and the inverse power of the length dependence for the signal), lending additional support to the astral stimulation model. Another, later model (Yoshigaki, 2003) added the central spindle hypothesis and tested all three mitotic apparatus–dependent models in the unified DE-based framework.

Subsequent experimental studies, inspired in no small part by these modeling efforts, led to the discovery of molecular pathways causing the activation of cortex contractility. These pathways include, but are not limited to, centralspindlin, a complex consisting of a kinesin-type molecular motor delivering the signal to the MT plus ends, and a Rho-family GTPase-activating protein triggering reactions resulting in actomyosin assembly and contraction (for a comprehensive review, see D’Avino ; Mishima, 2016). It would seem that the discovery of centralspindlin supports only the cortex activation models, that is, the astral stimulation and central spindle models. However, some studies, most prominently that of Canman ), indicated that MTs have simultaneously opposite effects: stable MTs boost, but dynamically unstable MTs suppress, cortical activation. How could that be? Again, one of the ways to answer this question is by the use of modeling.

This question was addressed by Odell and Foe (2008) by means of AB modeling. They hypothesized that there are two effects. 1) There is an indirect cortical relaxation effect, because centralspindlin binds to and walks along all MTs but rarely reaches the cortex on dynamically unstable MTs because they continually depolymerize before centralspindlin arrives. Dynamically unstable MTs thus sequester centralspindlin away from the cortex, suppressing actomyosin activation globally. 2) Some astral MTs, aimed primarily toward the cell equator, stabilize during anaphase and thereby become effective rails along which centralspindlin reaches the equatorial cortex. To simulate this scenario, a few hundred thousand agents were used: hundreds of MTs, each consisting of hundreds of polymer segments, thousands of centralspindlin complexes drifting on MTs and diffusing in the cytoplasm, and so on—even yolk particles to account for the excluded volume of the cytoplasm. Tens of microscopic processes were simulated, including transport and diffusion of centralspindlin, repeated growth and shrinking of MTs, nucleotide exchange on tubulin subunits governing MT dynamic instability, elastic bending of MTs, and much more. Massive simulations taking tens of hours on a powerful computer cluster resulted in stunning snapshots (Figure 2C).

Of interest, the authors failed to find a reasonable combination of parameter values for which these two effects caused a significant buildup of centralspindlin at the equator. This is because diffusion scattered centralspindlin as soon as the motors walked off the MT tips. Therefore the authors posited a third hypothesis: upon reaching the MT plus end, centralspindlin does not walk off but stalls at the tip. Simulation of the altered model achieved the desired result. The complex behavior captured by this model brings to mind a colorful analogy made by the grandfather of cytokinesis research, Ray Rappaport, who said that cytokinesis is not a beautifully made, finely adjusted Swiss watch, but rather an old Maine fishing boat engine: overbuilt, inefficient, yet never failing.

Before we start to compare DE and AB modeling in general, let us mention briefly that after we understand better how the spindle positions the furrow, the next question is how the spindle itself is positioned. Many studies, including modeling ones, have been devoted to this question, culminating recently in a combination of microscopy, micromanipulation, and AB modeling approaches that resulted in a model according to which dynein motors distributed both on the cell cortex and cytoplasmic scaffold pull on the astral MTs, creating a force balance that positions the spindle (Minc ). In addition, the fourth model, cortex self-polarization, has recently attracted quantitative modelers. Using reaction-drift-diffusion PDE systems to describe the dynamic densities of cytoskeletal and PAR proteins, they have found combinations of nonlinearities in the reaction terms and key feedbacks that could support self-polarization (Dawes and Munro, 2011; Goehring ), conceptually similar to the Turing mechanism.

Finally, there are fundamental questions about the assembly and mechanics of the contractile actomyosin cortex; in this case, too, modeling offers insight and possible answers. One of the earliest of these questions was asked by White and Borisy (1983): if there is an initial uniform distribution of elastic contractile elements around the cortex of the cell with high hydrostatic pressure inside and then the contractility is relaxed at the poles, would the model predict the observed sequence of cytokinetic cell shapes from one spherical to two daughter cells? A DE model gave a negative answer because the saddle-like curvature of the deepening furrow at the cell equator attenuates the effective force of scission. Note that such a DE model requires simulating mechanics equations, the physics of which is based on force balances (Danuser ). Mathematical forms of both force balance equations and reaction-diffusion-drift equations that we discussed earlier are very similar. Note also that other examples of AB modeling are the very useful (and very computationally expensive) molecular dynamics simulations based on solving mechanics equations describing molecules interacting with each other through realistic physical forces in aqueous media.

White and Borisy (1983) then realized that full cell division can be explained if contractile elements are pulled into the cell equator by the actin flow and align along the equator. Later, much more sophisticated DE models were able to mimic the mechanics of cell division with even greater precision (He and Dembo, 1997; Poirier ). Recently a detailed DE model of the actin cortex that includes signaling dynamics of Rho GTPases governing the cortical mechanics predicted highly nontrivial oscillatory excitable behavior at the onset of cytokinesis, which was confirmed experimentally (Bement ). Not only DE models proved to be useful to describe the mechanics of cytokinesis: an AB model (Vavylonis ) answered another fundamental question—how does a contractile structure self-assemble from a disordered array of actin and myosin filaments? This question asks for an AB modeling approach because the nature of this self-assembly is microscopic and random, involving a relatively small number of discrete agents characterized not only by positions but also by orientations and lengths. The emergent pattern formation, based on the coalescence of random search-and-capture–driven actin–myosin contractile units, is too complex to intuit and too multidimensional, heterogeneous, and noisy to be described by PDEs.

HOW TO CHOOSE BETWEEN DE AND AB MODELS

Consideration of two modeling studies described earlier (Devore ; Odell and Foe, 2008), both of them unquestionable but not absolute modeling successes, illustrates how to choose between two model types. Devore ) wanted to explore a conceptual rather than detailed model because molecular pathways of stimulating the furrow were not known at the time. Therefore the number of molecular species was minimal, explored geometries and dynamics were simple, and the resulting PDE were easily solvable numerically even at the time when computer simulations were much less routine. Perhaps the greatest advantage of the resulting model was the ease of scanning the parameter space. This said, their DE model was based on the unexamined approximation—that discreteness and randomness of the finite number of MTs, which are significantly splayed apart near the cell cortex, do not affect the continuous deterministic result. In addition, another assumption, that the signal delivered by MTs decreases as MT inverse length in some power, was relatively drastic and not supported by explicit modeling. If this study was done now, it would be a good idea to complement the DE model with the AB model to see how good the approximations that stem from these two assumptions are. The most desirable practice is to combine a more conceptual DE model with more detailed AB model, although this is not always possible, not to mention very laborious.

Odell and Foe (2008) made the decision to use AB modeling mainly because MT dynamic instability was hypothesized to be of primary importance for the astral stimulation mechanism and they wanted to develop a detailed model with a great number of processes with known or hypothesized rates included and they wanted to exclude approximations of truly random processes with deterministic ones. However, assuming a great number of MTs, one could easily investigate a DE model with MT dynamic instability accounted for by analytical expression for average MT density. Comparing predictions of such a simplified DE model with the full AB model would be very informative. Ultimately, the most important aspect one has to think about is what not to include into the model. For example, in the study of Odell and Foe (2008), one cannot help but wonder whether simulating MT elasticity and steric interactions, which consumed most of the computational time, was really necessary. Without those features, the model’s parameter space probably could have been explored better. These two studies illustrate that the choice between the modeling types is nontrivial and not unambiguous.

SUMMARY: ADVANTAGES AND SHORTCOMINGS OF DE AND AB MODELS

AB models are gaining popularity and for a number of reasons will likely become the dominant modeling tool in the future: one reason is that AB models are often considered to be in silico reconstitutions of biological systems because, if sophisticated enough, such models produce a life-like simulation of the cellular subsystem. The AB model allows computer experiments with in silico systems with an exquisite control of all parameters and ease of simulating biochemically or genetically perturbed systems. Because biological systems do actually consist of a large number of agents whose simple interaction rules produce mind-bogglingly complex behavior, the AB philosophy is very close to biology, and thus, often fewer approximations have to be made to build an AB model. An AB model is often much simpler than a DE model, especially in cases in which complex geometries, a high number of dimensions (including, besides spatial dimensions, distributions in angle and size), heterogeneity and anisotropy, and a great number of types of agents are involved. In addition, AB models capture stochastic effects more naturally than stochastic DE models. When a biological system consists of few discrete objects, the continuous approximation required for a DE model is not faithful. Last but not least, computer-savvy and quantitative-minded biologists can do AB modeling without the need to study mathematics.

However, of course, a number of catches ensure that DE models will never cease to be useful. The main problem with AB models is that it is often much harder than with DE models to get a qualitative insight: in silico systems tend to become too complex; in a way, we get interesting results but can only guess how the assumptions led to them. Fewer methods and software exist for AB models. Some artifacts are inherent to some types of AB models; for example, there are artificial oscillatory solutions in Boolean models and difference equations that do not correspond to reality. There are no analytical solutions to AB models, no benchmark cases. Usually, multiple simulations and vast and nontrivial statistics are necessary to extract meaningful insight from an AB model. It is usually hard to explore parameter space with AB models. Last but not least, AB modeling is hard to teach; unlike DEs, to which many textbooks are devoted, descriptions of AB modeling tools are scattered and not systematized.

All these minuses do not negate the pluses of thoughtful AB modeling. Many questions have to be considered before embarking on the hard modeling journey: Is modeling needed for this study at all (what is the question to be answered by it)? Do we honestly list and examine all assumptions and simplifications? Do we really look for a prediction, even if it contradicts the data, or are we subconsciously trying to confirm our preconceptions? Only when these questions are pondered can we start deciding on the method and type of the model.

COMPARISON OF DE AND AB MODELING (FIGURE 1)

I. Biology: Both AB and DE modeling start with formulating a hypotheses about the key players involved in the process in question. Each agent (e.g., protein, molecule, organelle, cell) interacts with other agents and/or the environment. Typical examples are biochemical (signaling, binding/unbinding) and physical (e.g., elastic collisions or excluded-volume effects) interactions. These interactions can change the position (transport effects of diffusion and directional drift) or the nature of the agent (chemical reactions). The model more often than not includes stochastic effects, reflecting variabilities in size, speed, or other properties of the agents, and random movements.

II. Model: A typical AB model consists of a large system of equations simultaneously updating positions and chemical states of each agent. Here (Figure 1, AB-Model), we show one such equation (Langevin equation in this case, but many other types of equations, or even equation-less particle rules, are often used), which updates the position of the ith agent in small time increments due to its stochastic Brownian motion (first line) and due to deterministic directional shifts as the result of interactions with all other agents (second line). For simplicity, we do not show reaction terms here. To translate the biological problem into a DE model, a continuous function describing the density of the agents has to be introduced, and reaction, drift, and diffusion terms containing partial derivatives of the density have to be derived. Many explicit and hidden approximations go into such derivations. In the model shown here (Figure 1, DE-Model), we exhibit the diffusion term and the term responsible for the drift and the reactions.

III. Results: AB models can be explored by performing many simulations using either available software or custom-made coding. Typical approaches include Monte Carlo or Gillespie algorithm simulation, in which many (often costly) simulation rounds are necessary to determine the mean behavior of the system. An advantage, however, is that one also gains knowledge about the variability of certain quantities (e.g., as depicted in Figure 1, AB-Model Results, the mean density of an agent).

Each problem at hand comes with a set of parameters (such as particle speed, reaction rates), and the system’s overall behavior might depend crucially on them. In Figure 1, we show aggregation and uniform distribution as examples of possible behaviors. For AB models, it is usually impossible to explore the whole parameter space (not only because of the computational cost), and therefore there is always a risk that a set in parameter space and the corresponding behavior are missed (compare AB and DE results in Figure 1). In other words, AB models rarely allow a full understanding of the system.

For DE models, on the other hand, the existing toolboxes include both analytical and numerical methods. In rare cases, analytical solutions for mean values can be obtained. Sometime a complete phase diagram—dependence of all possible patterns and behaviors of the system on the whole set of parameters—can be deduced (Figure 1, DE-Model Results). All this, however, can be done only for relatively simple DE models relying on assumptions that are sometimes hard to verify and on gross approximations.

IV. Note: Even though many of the mathematical models used in biology are AB or DE based, there of course exist combinations of those two approaches, as well as models employing other techniques. Examples include models using recurrence relations, graph theory, or integral and integrodifferential equations. In connection with models for swarming and collective motion, the Vicsek model (Vicsek ), a minimalistic AB model, has been used as a starting point for derivations of macroscopic, DE models (e.g., Degond and Motsch, 2008). Systematically deriving such continuous models from their microscopic counterparts is very meaningful because it allows harvesting the advantages of both approaches.