Setting the stage
A long-term biological goal for this project is to develop scientifically useful, validated simulations of leukocyte recruitment under physiological and pathophysiological conditions. The envisioned LEUKOCYTES, functioning within an analogue of the wet-lab context, will have many behaviors, properties, and characteristics that mimic those of referent in vitro and in vivo systems. Stated differently, there will be similarities between leukocyte phenotype and the "phenotype" of simulated leukocytes, as illustrated in Figure 1. The expectation is that increasing phenotype similarity will require, and can be achieved in part through, similarities in design plan and in generative mechanisms. Once they are achieved, the model will have become a scientifically useful analogue.
Our approach to designing and building the envisioned analogues is motivated by aspect-oriented software development and the middle-out approach first suggested by Brenner [22] and later detailed by Noble [23] for building models of complex biological systems. We adapted and merged these ideas in arriving at our working definition of middle-out modeling. First, specify the biological features under study (e.g., cell rolling and attachment to a surface). Each feature becomes an aspect of a software device – an archetype model. That archetype is iteratively transformed into a functioning analogue of the referent in vitro biological system. It is then validated and improved iteratively.
Following that approach, our first task was to pick a set of properties, p1, p2, ... pj (within a hyperspace of mechanistic properties) – around leukocyte rolling and adhesion in vitro. We specified Set A, Table 1. We next asked what software device could we implement to realize p1? How can we realize p2, etc.? We then asked, how can we realize p1- pj, all at the same time, using the same software device? Getting answers to those questions required exploratory modeling and simulation, occasionally using different modeling and simulation support packages. The in silico properties and the device we created to realize p1-pj formed a foundational analogue, the above archetype model, from which we intended to expand outward (up, down, and even sideways) to establish a reductive hierarchy of components that could be used to identify important systems biology principles.
Creating an analogue phenotype
The ISWBC is an abstract model of a leukocyte, and when placed in an appropriate, simulated environment, it is capable of mimicking a set of targeted characteristics. The ISWBC exists within a system that is a model of an in vitro flow chamber system (cells, media, the device in Figure 2B, atmosphere, etc.). The ISWBC is a multilevel composite model. The components are independent models of their referent leukocyte or environment counterpart. For example, the in silico component representing the MEMBRANE UNITS in a portion of the trailing edge of a ROLLING LEUKOCYTE is itself a model of the corresponding leukocyte feature.
The envisioned relationship between measured LEUKOCYTE behaviors and corresponding in vitro measures is illustrated in Figure 1. Pictured are overlapping (but not intersecting) sets. One set (Figure 1, large circle) contains the results of experiments that measured specific phenotypic attributes, properties, or characteristics of leukocytes in vitro (e.g., all three sets in Table 1). The smaller sets contain the results of simulation experiments that measure corresponding in silico phenotypic attributes, properties, or characteristics (such as Set A, Table 1). Our motivating hypothesis has been that, when attention is focused on the same features, and the measured behaviors of the two sets are similar, then there may also be useful similarities in the generative mechanisms of the two systems. Those similarities can be explored by iterative testing and refinement of the analogue coupled with related wet-lab experiments. An example: distance traveled by a LEUKOCYTE in silico (or a leukocyte in vitro) under various simulated (or actual) shear stresses. We cannot expect these measurements to overlap completely or even have the same units. Observations on behavior similarity can be used simultaneously to eliminate invalid model features and identify gaps in our knowledge of the experimental in vitro system.
The first step in transitioning from an archetype to an improved LEUKOCYTE was to thoroughly document the in silico properties and characteristics of the former. By doing so, we improved insight into its strengths and weaknesses. The next step was to obtain increasing overlap with measures of in vitro attributes. That was done systematically by iteratively revising together the hypothesis and the analogue by following the five steps below. An example hypothesis: the current analogue can acceptably simulate the first two attributes [p1, p2,] in Set A, Table 1. Figures 7 and 8 are examples of successful simulations. We begin the next iterative cycle with step 1.
1) Select an additional in vitro property or characteristic that is related to those already in the target set, and for which wet-lab experimental observations are available, such as the third property in Table 1. Identify the additional property as p3.
2) Add p3 to the targeted set, yielding an expanded set [p1, p2, p3], as illustrated by 2 in Figure 1B.
3) Determine if addition of the new attribute invalidates the current analogue, and if so, why. If not, repeat step 2.
4) Revise the model iteratively, possibly by adding mechanistic detail, until the measured phenotypic attributes of the revised analogue are sufficiently similar to [p1, p2, p3].
5) The next iterative cycle repeats steps 1–4 with p4, etc.
Our approach to achieving an envisioned ISWBC was strongly influenced and constrained by the ten capabilities listed in the Introduction. We reasoned that to achieve our long-term goal it would be essential for future ISWBC descendents to exhibit those capabilities.
Design
The conceptual model of the referent system, which is fundamental to the design of the in silico model, is as follows. Functionalities within the leukocyte membrane are grouped into a large number of similarly capable modular, functional units. Leukocyte rolling in vitro involves a balance between the hydrodynamic force exerted by fluid flow and the resistant forces resulting in part from transient bonds being formed between the leukocyte and the substrate. Shear force coupled with random events causes bonds to break at the rear of the contact zone. Bonds formed within the contact zone are by nature transient and will dissociate even in the absence of a tensile force. Only a few bonds need to be present at any time in order to maintain a rolling interaction. In order to maintain that interaction, broken bonds must be replaced by new ones formed elsewhere within the contact zone.
The transition from rolling to adhesion is mediated by the integrin receptors that, when induced into a high affinity state, are capable of forming effective ligand interactions that are sufficiently stable to withstand the hydrodynamic fluid force. However, they exist natively in a low affinity state and can only aid in supporting rolling interactions. There are many hypothesized mechanisms of integrin activation. We explored one of the simplest: a chemokine receptor, upon binding to its chemokine ligand, induces a local VLA-4 integrin molecule into a high affinity state. The resulting high affinity integrins form stronger and longer-lasting bonds with VCAM-1. Diffusion and clustering of VLA-4 integrins is also hypothesized to play a role in firm adhesion. However, to keep things relatively simple at this early stage, we have not simulated those events.
Achievements
The ISWBC components have been verified, plugged together, and operated in ways that may represent (at a relatively high level) mechanisms and processes that influence leukocyte rolling and adhesion. The approach has allowed us to represent the spatial and discrete event phenomena that are thought to occur during leukocyte rolling, activation, and adhesion. It has allowed us to represent apparently stochastic elements within the system that may play a vital role in determining leukocyte behavior in vivo and in vitro. Our earliest archetype had a CONTACT ZONE comprised of 2 × 3 MEMBRANE UNITS (48 for the total MEMBRANE of the LEUKOCYTE). With it, we failed to even approach the first two targeted attributes in Set A, Table 1; we were unable to simulate the stochastic fine structure present in the in vitro data because the discretization was too coarse. It is common to observe such discretization artifacts during simulation development when the granularity is too coarse. To better achieve the targeted behaviors, we decreased the relative MEMBRANE UNIT size (decreased granularity) to represent the CONTACT ZONE by 3 × 6 units and then in steps to the current size of 8 × 10 units (600 UNITS for the total MEMBRANE). With each increase in granularity, simulated results more closely mimicked the in vitro data, as judged by inspection. We also explored even finer grained CONTACT ZONES, up to 20 × 30, without further gains in similarity. Taken together, these results suggest that membrane functionality within and adjacent to the contact zone is distributed into many (from our results, about 80) equally capable, somewhat autonomous units.
Though the ISWBCs are relatively simple, we have demonstrated that they can generate the targeted attributes under three different experimental conditions. Traditional models typically focus on just one condition. ISWBCs mimic the dynamics of leukocytes rolling separately on a P-selectin or VCAM-1 substrate. In addition, they mimic the transition from rolling to adhesion on P-selectin and VCAM-1 in the presence of GRO-α chemokine. Importantly, when simulating populations of leukocytes under different experimental conditions (combinations of substrate molecules), the system generates quantitative population-level data that are similar to data from in vitro experiments.
One hypothesis of leukocyte activation suggests that leukocytes roll along the vessel wall progressively engaging with chemokine signals to reach a global activation threshold. Once achieved, fully activated integrins can then support firm adhesion. However, there is growing support for an alternative hypothesis: leukocyte adhesion induced by immobilized chemokines occurs via local and spatially restricted signaling to nearby integrins. Our results support the latter hypothesis. In an elegant study, Shamri et al. [24] showed that leukocyte engagement with chemokine locally and reversibly induced the β2 integrin, LFA-1, into an intermediate affinity state that was essential for interaction with ICAM-1 under flow conditions. Concurrently, engagement with ICAM-1 induced the signaling events needed to fully activate local LFA-1 to the high affinity state needed to support firm adhesion. Because of the modeling approach used (providing the ten capabilities, etc.), it will be relatively easy for future descendents of the ISWBCs to represent such fine-grained, cooperative behaviors and spatially restricted events.
Appraisal of the model assumptions
We have simulated human leukocytes of slightly different sizes (neutrophils, lymphocytes, monocytes), however we kept the CONTACT ZONE the same for all simulations. It has been observed that leukocytes of greater size have a greater area of surface contact, and thus can form more bonds that are adhesive. However, biomechanical analysis shows that with increasing leukocyte size there is an increase in the amount of force and torque experienced. This higher magnitude is believed to translate to higher forces experienced by bonds. Cozen-Roberts [25] presented an analysis that was verified in vitro, which predicted that the disruptive force increases faster than the pro-adhesive effect of increased contact area. Rather than changing the dimensions of the contact area when simulating experiments with different leukocyte types and sizes, we fixed that dimension and compensated by using greater RearForce values to represent the same shear force when simulating larger leukocytes (see Table 5).
The conceptual model above is widely accepted for neutrophils. However, it is not entirely known how well it applies to other leukocyte types, such as lymphocytes and monocytes. Varying receptor types and densities found on the different leukocyte types may be a factor contributing to differing mechanisms of adhesion. For example, lymphocytes and monocytes express basal levels of VLA-4 integrins sufficient to support rolling [17]. In contrast, it has been shown that transmigration across endothelium or exposure to chemotactic stimuli is necessary to induce expression of VLA-4 integrins on human neutrophils [26]. Therefore, rolling and adhesion through the VLA-4 integrin may require additional or different mechanisms for neutrophils than for lymphocytes and monocytes. Our objective during this project was not to construct a comprehensive model of rolling, activation, and adhesion for all leukocyte types. However, when simulating monocyte rolling and adhesion, we assumed that monocytes and neutrophils share similar rolling characteristics on P-selectin. Likewise, we assumed that the characteristics of monocyte rolling on VCAM-1 are similar to those for lymphocytes. The only difference was presumed to be the densities of each receptor type found on each leukocyte type. Parameter values representing the densities of the VLA-4 integrins for lymphocytes in experimental condition 2 and monocytes in experimental condition 3 were different to reflect the different VLA-4 density values reported in the literature. Values for PSGL-1 sites/human monocyte were not found in the literature and were assumed to be similar to values reported for human neutrophils.
Comparison to other leukocyte rolling and adhesion models
The discrete-time models used in the Adhesive Dynamics simulations by Hammer and co-workers are the most developed models of leukocyte rolling and adhesion to date [9–11]. Leukocytes in their models are idealized as solid spheres decorated with rod-like microvilli containing receptors at their tips. Their simulations have allowed them to explore the molecular properties of adhesion molecules, such as reaction rates and bond elasticity, and how these properties may relate to macroscopic behavior such as rolling and adhesion [9, 11].
In their simulations, the position of the cell is determined at each time step from the net force and torque on the cell using a hydrodynamic mobility function for a sphere near a plane wall in a viscous fluid. The net force and torque acting on the cell from bonds, fluid shear, steric repulsion, and gravity are calculated. When calculating the forces on the sphere due to fluid shear, their model uses the Goldman equation [27]: the amount of force a sphere experiences from an applied shear force is calculated by assuming that the sphere is solid. At each time step in the Adhesive Dynamics simulation, positions of bonds on the spherical particle are tracked enabling the authors to calculate the forces that each bond experiences.
It was not our intention to discover ISWBC parameterizations that would yield simulation results that tightly fit specific sets of experimental data. Our approach has been focused on simulating a targeted set of system-level properties (Figure 1). Therefore, we did not, for example, attempt an explicit mapping between the shear force and the amount of force experienced by bonds at the rear of the leukocyte. Leukocytes are notoriously deformable and therefore we elected to avoid a mapping such as the Goldman equation. In addition, the mapping from the force on the cell to the force on the rear bonds can be quite complex, as PSGL-1 and VLA-4 are both ligands that are concentrated at the tips of stretchy and heterogeneous microvilli. It should be noted that shear stress for selectin tethers at the rear of neutrophils have been estimated previously using the Goldman equation to be 124 ± 26 pN per dyn/cm2 [28].
To account for the effect of an applied force on the kinetics of bond dissociation, Hammer and co-workers have employed both the Bell model and Dembo Hookean Spring model (see [2]). The Bell Model predicts bond dissociation rate as a function of applied force, whereas the Dembo model treats bonds as Hookean springs and relates dissociation rate to the length of the stretched bond. Use of a similar fine-grained representation of bond dissociation may have yielded more precise simulation results, but that level of resolution and precision was not needed to meet the simulation objectives in Table 1. Rather, we used a simple model to relate the probability of bond dissociation with bond force (Figure 4A); it was motivated by in vitro experiments by Park et al. [15], and Zhang, et al. [29]. Because the time resolution within the ISWBC is low (a tenth of a second), using the linear approximations in Figure 4A proved reasonable: at small bond force values, we assumed a linear relationship between the probability of bond dissociation and bond force. At larger bond force values, we reasoned that the rates of bond dissociation for ligands are too fast to be distinguished on our time-scale. Consequently, BONDS experiencing larger forces are represented as having undergone a DISSOCIATION event during a simulation cycle with a probability of 1.
Many molecular level details are believed to impact effective bond formation and breakage within small portions of a leukocyte membrane and the surface. Examples include contact irregularities, local dynamics, including ligand relocation within the membrane, force history of bond loading, and bond compliance. All are aggregated and controlled in the ISWBC by an event probability. When explanation of system level behaviors requires a more detailed representation, one or more of these factors can be specifically represented, without compromising the function of other ISWBC components. The probability parameters will remain, but their values and explanation will have changed.
In the most recent version of their model, Hammer and co-workers represent the p38 mitogen-activated protein kinase (MAPK) global integrin activation mechanism leading to leukocyte arrest induced by E-selectin engagement with its ligand PSGL-1 [30]. A simple modular Hill function is used to represent the overall MAPK cascade leading to integrin activation. At the end of each time step, the total number of E-selectin-PSGL-1 bonds are counted and used as an input to this function. We have chosen to explore an alternative integrin activation mechanism in which detection of chemokines activates local integrins. Future embodiments of this modeling effort may include the E-selectin-dependent signaling pathway.
Recent experiments have shown that leukocyte rolling on substrate-coated, flow chamber surfaces is significantly slower and smoother than that of microspheres coated with ligands [15, 31]. It has been hypothesized that the deformability of leukocytes can influence rolling by increasing the contact area. In vivo, leukocytes rolling at blood wall shear rates of 800 s-1 elongate in the direction of flow to 140% of their estimated undeformed diameter and have a 3.6-fold increased contact area with the endothelium than at blood wall shear rates of 50 s-1 [32]. These experiments have influenced other modelers to explore the mechanical and rheological properties that influence leukocyte morphology and deformation. The 2-D elastic ring model [33] and the 2-D compound drop [34] model offer an account for the deformability of leukocytes during rolling. However, because of the way receptor-ligand interactions are treated deterministically, these models do not produce the characteristic jerky stop-and-go behavior of rolling leukocytes.
Recently, Jadhav et al. created a 3-D model that provides a representation of leukocyte deformability: it can mimic the stochastic nature of receptor-ligand interactions [35]. The model uses a mathematical formulation called the Immersed Boundary Method to simulate the measured motion of a leukocyte by representing it as an elastic capsule in a linear shear field. They use a Monte Carlo method for simulating receptor-ligand interactions, enabling them to simulate the jerky stop-and-go movement of rolling leukocytes. Their model explored the effects of membrane stiffness on cell deformation and leukocyte rolling. It is believed that the contact area between the surface and the leukocyte membrane would increase for leukocytes under increasing shear rates. By introducing deformation in their model, they are able to observe this change in contact area. Our model currently ignores this phenomenon: we kept the contact area fixed for all simulations. Future ISWBC descendents can, when needed, allow each CELL to vary its own contact area.
Considerable knowledge has been gleaned from these different models of leukocyte rolling and adhesion. They have provided insight into how molecular parameters of adhesion molecules or cell deformability may influence leukocyte rolling and adhesion. The knowledge gained reinforces the merit of investigating complex phenomena using different modeling methods. Our approach is quite different and is aimed at studying the interactions of the components and possible mechanisms that may give rise to emergent leukocyte- and systems-level phenotypes.
Computational models of cell motility
Beyond the primarily mechanical models considered above, what other modeling methods might one consider to achieve the goal presented above under Setting the Stage? Large-Q Potts models are based on extensions of cellular automaton (CA) models [36], including lattice-gas CA, with features adapted from the Ising model. They have been used, for example, to study the cell sorting of embryonic cells [37], and extended to simulate a variety of different biological phenomena [38]. A further extension using a layer for partial differential equations (describing cAMP dynamics) enabled elegant simulations of the complex cell movements occurring during culmination of the morphogenesis of the slime mould Dictyostelium discoideum [39]. Importantly, the simulated cells were represented as a group of connected automata making the basic model scale subcellular. Important details of these methods are collectively reviewed by Alber et. al [40].
Such extended Potts models have been the preferred approach for describing cell motility at a subcellular level, including aspects of cell shape, surface chemistry, and internal structure. Nevertheless, Meyer-Hermann and Maini found it necessary to develop an alternative, agent-oriented model architecture to interpret recent the two-photon imaging data of lymphocyte movement within lymph nodes [41]; the jerky and individualistic nature of the imaging data is similar to that of leukocyte rolling. They refer to their system as the hyphasma (the Greek word for tissue) model.
The hyphasma LYMPHOCYTE is a connected collection of subunit objects arranged on a 2D square grid centered on a virtual barycenter. Each one has a polarity vector, a list of its constituent subunits, internal velocity states, and a CELL volume that determines the number of subunits. During movement, the barycenter is shifted in the direction of the polarity vector and border subunits are randomly moved to lattice points near the newly established barycenter. The simulated results exhibit Capability 4 (Turing Test) and support the hypothesis that lymphocyte movement within secondary lymphoid tissue may not be a consequence of chemotaxis or haptotaxis, but a simple random walk. It seems likely that the hyphasma model could be used to simulate data like that in Figures 5, 6, 7, 8. However, extending the model so that it delivers Capability 3 (Mapping) seems problematic. Relative to the abstract nature of the subunits, the authors state: "a force balance equation for the cell subunits is not necessarily a correct description of active processes of deformation and reshaping of cells. Thus, the interpretation of cell subunit velocities in terms of forces has to be considered as an approximation to more complex internal processes within the cell." Simulating data like that in Figures 9 and 10 is beyond hyphasma's current capabilities, but then the model was not developed with the intent that it be adaptable and extensible to represent such data (Capabilities 5 and 8). As stated in Introduction, when developing the ISWBC we had in mind objectives beyond simply simulating the in vitro data presented. We sought and presented a method for assembling individual components (that map to identifiable biological counterparts) according to a design, and then showed that the constructed ISWBC can exhibit behaviors that match those observed.
Future directions
The ISWBC described here is an important first step towards our long-term goal of building in silico devices that can simulate leukocyte behaviors in a variety of in vitro experimental systems even in the face of considerable uncertainty. They are useful objects for experimentation, as are the referent in vitro models. The current ISWBC is capable of mimicking only a few (Table 1) of a long list of desired phenotypic attributes of leukocytes observed in vivo and in vitro. The expectation is that ISWBCs can be iteratively refined to become increasingly realistic in terms of both components and behaviors. Because we strove to deliver the ten capabilities listed in the Introduction, any of the above, abstract, low resolution components can be replaced with validated, more realistic, higher resolution composite objects that map to more detailed biological counterparts, without having to reengineer the whole system, and without having to compromise already validated features and behaviors.