Modeling emergent tissue organization involving high-speed migrating cells in a flow equilibrium

There is increasing interest in the analysis of biological tissue, its organization and its dynamics with the help of mathematical models. In the ideal case emergent properties on the tissue scale can be derived from the cellular scale. However, this has been achieved in rare examples only, in particular, when involving high-speed migration of cells. One major difficulty is the lack of a suitable multiscale simulation platform, which embeds reaction-diffusion of soluble substances, fast cell migration and mechanics, and, being of great importance in several tissue types, cell flow homeostasis. In this paper a step into this direction is presented by developing an agent-based mathematical model specifically designed to incorporate these features with special emphasis on high speed cell migration. Cells are represented as elastic spheres migrating on a substrate in lattice-free space. Their movement is regulated and guided by chemoattractants that can be derived from the substrate. The diffusion of chemoattractants is considered to be slower than cell migration and, thus, to be far from equilibrium. Tissue homeostasis is not achieved by the balance of growth and death but by a flow equilibrium of cells migrating in and out of the tissue under consideration. In this sense the number and the distribution of the cells in the tissue is a result of the model and not part of the assumptions. For purpose of demonstration of the model properties and functioning, the model is applied to a prominent example of tissue in a cellular flow equilibrium, the secondary lymphoid tissue. The experimental data on cell speed distributions in these tissues can be reproduced using reasonable mechanical parameters for the simulated cell migration in dense tissue.


Introduction
There exists a number of agent-based off-lattice models to simulate tissue [1][2][3][4][5][6][7][8]. Many of them describe tissue organization and pattern formation [9] by the dynamics of growth and death of cells. Other 1 models are suitable for describing migration of a constant number of cells. The traditional research field of pattern formation involving migrating cells [3,4,10] is mostly based on a balance of proliferation and cell death. Cells are migrating with relatively slow speed with respect to other processes like diffusion of soluble substances. Approximations based on slow cell migration are justified even for such highly dynamic systems as the epidermis [8,11].
The model introduced here is designed for tissue in a flow equilibrium of high-speed migrating cells as found frequently in immunological tissue. Cells can enter and leave the tissue according to some dynamics imposed by the considered biological system. Feedback mechanisms between the cells and the substrate are incorporated because homeostasis in such systems is often established by feedback interaction [12][13][14]. This is modeled as transient excitation of a substrate or a small sessile cell population. The excitation is thought to couple back to the cell migration thus allowing for a regulatory system in a stable flow equilibrium.
The model is based on a regular triangulation [15,16] that provides the cell neighborship topology which is changing rapidly due to fast cell migration. The dual Voronoi tessellation [16,17] provides information about the shape of the cells including cell contact surfaces [8,11,18]. The physical properties are based on previously introduced short-ranged elastic interactions and actively generated intercellular forces [3-5, 8, 18]. Each cell is modeled individually and can incorporate a set of molecular properties or internal states that is appropriate for the system under consideration. Such properties may have impact on cell migration or cell-cell and cell-substrate interaction. The deterministic internal cell dynamics exhibit memory in the sense that, for example, internal variables representing the cell state like delay times are included. Cell migration under the influence of chemotaxis is also treated deterministic. However, cell orientation is stochastic for unguided random cell migration.
In order to illustrate the model features, the model is used to investigate certain features of cell behavior that are relevant for the formation of primary lymphoid follicles (PLF). These can be found in secondary lymphoid tissue (SLT) like the spleen and lymph nodes [12,14,[19][20][21][22]. Two lymphocyte cell populations of the PLF show fast migration: B cells and T cells. Two sessile populations, namely follicular dendritic cells (FDC) and fibroblastic reticular cells (FRC) serve as a complex 3D substrate for lymphocytes. The morphology of SLT is characterized by PLF and T zones. The PLF contains the B cells and the FDC. It is adjacent to the T zone harboring T cells and FRC. The functionality of the model features is demonstrated by analyzing the migration of the lymphocytes and their interaction with chemokines. Also the generation of a tissue in flow equilibrium is studied using a toy model inspired by PLF formation. On the basis of biologically proven interactions the model can explain a clear separation of T and B cells forming a T zone and an FDC containing follicle, as found in real tissue. However, more realistic models of PLF remain to be developed in future work.
The presented simulation technique has a direct connection to experimental constraints and exhibits the potential of building physically concise models of cell migration in tissue. It's novelty with respect to previous work relies on the focus on tissues composed of rapidly migrating cells which exhibit cell-flow-equilibrium. This extends the range of applications of agent-based off-lattice techniques to lymphoid tissue.

Method
The tissue formation problem of migrating cells is simulated using an agent-based model on top of a regular triangulation [15,16,23,24]. The regular triangulation is used to provide the neighborhood topology for the cells that allows for a continuous representation of cell positions and sizes in contrast to grid-based methods. The regular triangulation also provides information about the cell contacts and contact areas, and, within the limits of an approximation, about cell volume and shape [8,25]. The regular triangulation will be briefly summarized in the next subsection.
The simulation of cells is realized in a 3-level multiscale model. The first level is the internal state of the cells representing the dynamics of the phenotype of the cell (see Sec. 2.2). The second level models the contact interaction between cells including mechanical interactions with the environment and exchange of signals by membrane bound molecules (see Sec. 2.3). The third level incorporates long range interactions via diffusive substances for example chemotaxis, i.e. directed motion of cells upwards a concentration gradient which is induced by molecular chemoattractants (see Sec. 2.4). Thus, chemoattractants are derived from cellular sources and feed back to cell migration.

Regular triangulation
Each cell is represented as a sphere at position x with radius R. A vertex is defined by the pair X = (x, R). The regular triangulation is defined using the empty orthosphere criterion [15][16][17]. In three dimensions four vertices A, B, C, D forming a tetrahedron uniquely define an orthosphere (Fig. 1). The orthosphere is empty if for any other vertex V holds provided that the four vertices are oriented positively, i.e. if The symbols a x , a y , a z denote the coordinates of vertex A at the position a = (a x , a y , a z ) and R a the corresponding radius of the elastic sphere associated with the vertex. The notation for the vertices B, C, D, V is analogously defined. A set of non-overlapping tetrahedras covering a set of vertices forms a regular triangulation if all orthospheres attributed to these tetrahedras do not contain any further vertex. The regular triangulation, and with it the neighborhood topology, changes when cells are moving, and when they are added or removed from the system due to cell flux, cell death or cell proliferation. The corresponding algorithms have been developed and published previously for serial [23] and parallel computer architectures [24].

Internal cell dynamics
The phenotype of a cell is described by a set of internal cell variables φ. These include internal times to indicate when which type of event may happen. An example is the persistence time T p of cells during chemotactic motion [26][27][28][29][30]. When the persistence time T p has past, the cell can reorient to the local chemoattractant field [31,32]. Upon contact two cells can exchange signals via the contact surface. In a simple approach with unpolarized ligand/receptor distributions the signal strength is proportional to the contact area. The contact area is computed as the minimum between a sphere overlap and the common Voronoi face of these cells (Fig. 2). This choice relies on the fact that the Voronoi-contact area is the better description for high density of cells while the virtual overlap of the spheres is more realistic for low density systems [8,18]. In both cases the alternative measure for the shape leads to larger estimates which justifies the present choice. In general the minimum of the two contact areas is close to the realistic contact of two elastically deforming spherical cells.
Finally, a set of variables describes the mechanics of the cell: velocity, orientation, cell volume, and elasticity. These variable couple directly to the next level of description, the contact interaction of cells, and have no direct influence on other internal states but may influence the corresponding variables of neighboring cells.

Equations of motion
The contact interaction of cells is predominantly given by mechanical interactions. It is described by Newtonian equations of motion in the overdamped approximation [2,4,8]. In this approximation acceleration of cells and consequently conservation of moment can be ignored. With The active forces F act on cell i (if any) at position x i depend on the internal state φ i of the cell and the internal state of neighboring cells φ j , while the passive forces F pass depend on the positions of the cell i and its neighbor cells j ∈ N c i in contact with cell i. Active forces can also be exerted directly to the surrounding medium. Therefore the set N i includes all Delaunay-neighbors j independent of whether they are in physical contact with the cell i or not. All these forces are counter-balanced by the velocity-dependent drag forces F drag . Because of the overdamped approximation this results in an ODE system of first order for the cell positions.
The passive forces are composed of forces stemming from the cell's elasticity and compressibility where F JKR ij recollects elastic interactions and surface energy.

Elastic cell-cell interaction
Elastic forces between cells are treated according to the JKR-model [33]. They depend on the virtual with the cell's elastic moduli E i and E j , Poisson numbers ν i and ν j , and the surface energy σ ij . The force acts in direction of the normalê ij on the contact face The applicability and the limits of the JKR-model for modeling elastic interactions between spherical cells have been quantified recently [34].

Many-body interactions
The JKR-model defines a two cell interaction which exhibits an equilibrium distance mediated by the balance of surface energy and elastic repulsion. For two cells this leads to a negligible deviation of the volume attributed to the vertex and the relaxed volume of the real cell. However, cells will frequently interact with several cells such that specific situations can occur in which cells are strongly compressed by surrounding cells without a corresponding large relaxing force being generated. Thus, cells might remain in a highly compressed state for too long times because of the neglect of manybody interactions. To account for the cell volume and approximately ensure volume conservation, a cell pressure concept is included [4,8]. The pressure of cell i is calculated as deviation of the actual cell volume V i from the target volume V * where a linear compression model with compressibility K i is used. The forces resulting from this pressure are exerted between cells by adding the term to the passive cell forces F pass ij (x i , x j ) in (4). a ij is the contact surface of the cells with neighbor cells or the medium. In case of contact surface with the medium a ij is calculated as minimum of the Voronoi face and the sphere overlap with an equal-sized virtual cell being at equilibrium distance with respect to the JKR forces. To get the actual volume the minimum of the sphere volume and the Voronoi cell volume is taken which is motivated by similar arguments as for the contact area of two cells [8,18].

Force generation of migrating cells
Previous models that incorporate the active forces of migrating cells [2][3][4] use a net propulsion force that is exerted into the direction of motion by pulling on neighbor cells. However, cells also exert forces perpendicular to the direction of motion which can be even larger than the net force as it has been shown for keratocytes [35,36]. The target application is lymphoid tissue which requires to describe lymphocytes. However, the mode of migration is different for these cells. It depends on the so-called constriction ring [37][38][39][40][41]. Therefore the active force generation model of fast migrating cells is based on the constriction ring model thus taking into account the lateral forces during active migration of lymphocytes.
Cells use a ring to attach a part of their membrane and cytoskeleton to the extra-cellular matrix and exert outward directed pressure to their environment by contracting the rear of the cell (Fig. 3). The ring remains fixed with respect to the extracellular matrix and therefor moves towards the end of a cell during migration. In a simplified approach a new ring is generated at the front of a cell only when the preceding ring has reached the end of the cell. The force acting on cell i by exerting active forces on a neighbor cell j reads with the cell orientation o i , cell surface contact point x * ij , constriction ring center x * i (which generally differs from the cell's center), interaction area a ij and the pressure p * i actively exerted by cell i. Note that p * i is a parameter of the model and shall not be confused with the pressure in (7).
x i x ij o i * * Figure  The motion of cells can also be characterized using a motility coefficient that is similar to diffusion coefficient in Brownian motion. However, it is not just a result of physical contacts with the environment but also involves intrinsic features of cell migration. The motility coefficient D is influenced by the cell's speed v and the mean free path length [42]. In cells the latter one is determined by the persistence time T p during which the cell maintains its polarized structured fixing the direction o i of force exertion and therefore the direction of motion [26][27][28][29][30].
Thus the persistence time is a parameter describing the migration of cells which is required in addition to the force generation.

Friction forces
The component of the friction forces between two cells i and j in contact read taking into account only the tangential parts of the relative velocity consistent with the treatment of the cell-cell interaction by the energy conserving JKR-model (5). The force depends on the relative cell velocityẋ ij =ẋ j −ẋ i . The overdamped approximation and low Reynolds-numbers of cell migration justify the linear velocity-dependence [2,4,8]. The unit vectorê ij is the normal vector of the plane of the cell contact. The friction coefficient γ ij has the dimension of a viscosity times a length scale. γ ij is chosen proportional to the contact area between the cells. The total drag force is the sum of all friction forcesF drag ij with neighbor cells plus a term of the free surface of the cell having interaction with the surrounding medium. The overall drag force is then given by with medium viscosity η med and the cell-specific viscosities η i . A i = j∈N c i a ij is the surface in contact with other cells. A tot i = j∈N i a ij is the total surface of a cell. The friction coefficient fulfills the symmetry γ ij = γ ji . The form of the friction coefficients is motivated by the Stokes relation for the friction of a sphere at velocity v in a medium with viscosity η The geometry related factor 6π is absorbed in the cell-specific viscosity η i in (13). The coefficients η med and η i are chosen such that reasonable active forces have to be generated by the cells in order to match the cell velocities measured experimentally (see Table 1).

Reaction-diffusion system of chemoattractants for long-range interactions
The chemotaxis of cells is described by coupling the direction of migration o i of a cell i to the local chemoattractant gradient. According to the observation that leukocytes tend to have a persistence time T p between subsequent orientation changes [20,[26][27][28][29] the gradient is sensed by the simulated cells periodically and o i is kept constant in between. The concentration of a chemoattractant is computed solving the time dependent diffusion equation. The time dependence is included for two reasons. First, to account for a dynamics of sources that frequently generates new or removes present sources. The time to equilibrate the concentration after changing the source is far longer than the typical times for cell migration such that the cells can sense this temporal concentration change. Second, the internalization dynamics presented in the next section acts on a time scale comparable to the cell migration but much faster than the time scale required for the chemokine diffusion to be equilibrated on the tissue scale.

Internalization model
The internalization of chemoattractant receptors occurs naturally when they are bound by their ligand [43][44][45]. Studies showed that desensitization of chemoattractant receptors exposed to high chemoattractant concentrations [46] occur probably by internalization. Alternative possibilities would be a chemical modification of the receptors which is not considered here. Of note, some experiments  [47]. Thus, the presence of receptors is necessary but not sufficient to cause chemotaxis of cells. Most likely the suppression of the receptor function is mediated by a desensitization mechanism either by internalization or cross-desensitization [44]. Desensitization of cells with respect to a chemoattractant c is achieved by internalization of the receptors which have bound the chemoattractant with rate k i [45]. Thus the receptor comes in three states: free on the cell membrane (R), bound to the ligand on the cell membrane (R b ), and internalized with the ligand (R * ) (Fig. 4), leading to the dynamics: The binding of the ligand is characterized by the on and off rate constants k on and k off leading to the ligand-receptor complex R b . The rate k r describes the recycling of the internalized receptor R * into the free membrane form R. The basic assumption is that the total receptor content R tot = R+R b +R * is conserved, i.e. no terms describing the transcription or degradation of the receptor are included in (15). This view is supported by the fact that internalized ligand-receptor complexes dissociate such that free receptor become available in the cell [45]. The assumption allows to eliminate the equation for R b from the system (15).
The reaction equations (15) are completed by an unspecific decay κc of the chemoattractant, a source term Q, and a diffusion term ∆c for the ligand. The decay term is in accordance with the fact that chemoattractants are inactivated and processed by all cells without reaction with the proper receptor [44]. This limits the life time of a chemoattractant molecule. The receptors are also transported by cells such that the whole system (15) is a reaction-diffusion system. It is solved using a splitting method of first order in time [48,49].

Model of primary lymphoid follicle formation
This section briefly describes the underlying biological concepts of PLF dynamics. The required ingredients to study the formation of the PLF are the B and T cell flow and the dynamics of the sessile cell populations, i.e the behavior of FDC and FRC. The parameters for the simulation are given in table 1. A detailed parameter estimate for the reaction-diffusion system (15) is given in the appendix.

Flux of B and T cells
B cells are constantly entering SLT predominantly via specialized blood vessels (reviewed in [21]). It is less clear how B and T cells leave this tissue. Recent experimental data suggests that lymphatic vessels guide cells to leave SLT [21]. In a simplified model cells enter and leave the tissue through separate spherical areas that do not further interact with the system. However, one should be aware that in real SLT the lymphatic endothelium which may be involved in the efflux of cells is located exclusively outside the PLF [50][51][52][53][54][55][56][57][58][59][60]. Thus, there may be some interaction that couples PLF formation with the location of exit spots which we ignore in the presented toy model. In order to allow the cells to remain and interact within the tissue for a certain time they are not permitted to use the exit areas before a minimal transit time has past. This time is about 3 hours according to a receptor dynamics [61]. This is consistent with measurements of the minimal transit time of lymphocytes [52,62,63]. In order to establish a sufficient cell efflux a chemoattraction towards the exit spots is assumed. This is supported by experimental evidence [21,61,64]. Due to a lack of information the corresponding receptor is treated with a simplified model without taking into account receptor internalization.

Migration of B and T cells
The migration of B and T cells within SLT follows similar principles. Both cells can perform chemotaxis and random motion. It has been shown experimentally that B cells are predominantly attracted by the chemoattractant CXCL13 [21]. Both, B cells and T cells respond to CCL21 with T cells being more sensitive [21]. The corresponding receptors are CXCR5 (for CXCL13) and CCR7 (for CCL21). FRC are the source for CCL21 and FDC produce CXCL13 (FIG. 5).

Dynamics of sessile cell populations
Although the origin of FDC is not clear there is strong evidence that they are derived from FRC found in the T zone of SLT [65]. Therefore, it is assumed that FDC are an excited state of FRC. The excitation is mediated by the contact signal LTα 1 β 2 provided by the B cells [12][13][14] (FIG. 5). The FDC-B cell interaction has been shown to be reversible [12] in the sense that the absence of the B cell signal lets the FDC vanish. In the model context that is described as the decay of the excited FDC state back into the FRC state. Note that, alternatively to the assumed scenario FRC and FDC may also have a common progenitor which can develop in either FRC or FDC depending on external stimuli. This is not considered explicitely in the present model. Within the PLF model an FRC differentiates to an FDC when the signal threshold for LTα 1 β 2 has been exceeded for a given time T FRC←FDC . The signal is determined by summing up all LTα 1 β 2 contributions from neighbor cells, i.e. surface density of LTα 1 β 2 times contact area. The differentiation is then instantly performed changing the internal cell states of an FRC into that of an FDC, i.e. replacing a source for CCL21 by a source for CXCL13. In a similar manner FDC differentiate back to FRC after the LTα 1 β 2 signal is below the threshold for a critical time T FDC←FRC . For simplicity and considering the lack of corresponding experimental data T FRC←FDC = T FDC←FRC and equal LTα 1 β 2 thresholds are assumed.

Sequence of events for follicle formation
The overall picture of PLF formation is as follows. A background of immobile FRC produces CCL21. B and T cells enter the tissue through the blood vessels which are represented by few (<10) spheres of 30 µm diameter [66] randomly scattered in an area of 150 µm size. When sufficiently large B cell aggregates form, FRC are induced to become FDC by LTα 1 β 2 thus replacing a production site for CCL21 by a CXCL13 source. This attracts more B cells and enlarges the forming PLF. T cells are kept outside the PLF just by coupling to the CCL21 which is produced around the PLF by the remaining FRC. Both, B and T cells, leave the tissue through exit spots which are represented as spheres acting as sinks for cells. Any dynamics on the exit spots that may result from PLF formation are ignored (compare Sec. 3.1).

Biomechanics of lymphoid tissue
The active forces exerted by lymphocytes to the ECM (F act i (φ i )) and to other cells (p * i , (10)) are not known explicitely (there are just estimates based on experiments of other cell types [35,[82][83][84]). The same applies to the friction parameters (η i and η med , (13)). For simplicity it is assumed that η i = η med for all lymphocytes. As the known average speed is essentially determined by the ratio of active forces and the friction parameters only 2 parameters per cell type remain unknown. These can be determined from the measured speed distributions for B and T cells within SLT [26,29]. A reduced cell system with a constant number of lymphocytes that respond to a chemotactic point source is used. The cell constituents are either B or T cells only or an even mix of both cell types.
The parameters for the JKR-forces are relatively well known. The parameters entering (5) are the elasticity of lymphocytes E i = 1kP a [34,[70][71][72], the Poisson number ν i = 0.4 [73,74], and the surface energy σ ij = 0 − −0.3nN µm −1 [75,76]. Thus the elastic cell interactions can serve as reference forces limiting the reasonable range of friction parameters. Therefore the fit to the speed distribution was performed with two qualitatively different friction regimes. In the low friction regime (η i = η med = 50 nN µm −2 s) the active forces that are required to produce the correct speed distributions are of the same order as the JKR forces. In this regime it turns out that the mode of measuring the cell speed becomes crucial: The speed distribution of the real system was measured on the basis of the average displacement over an interval of ∆T = 10 − 15 s [26,29]. in vivo [25] in vivo [28] low friction, ∆T=15s high friction, ∆T=15s  faster sampling with ∆T = 1 s the speed distribution becomes different and exhibits the influence of JKR forces that lead to fluctuations at rather high speed (FIG. 6 top right panel). The JKR forces dominantly lead to reversed or perpendicular motion of cells with respect to their migration axis when cells collide due to active movement. Thus the net-displacement by JKR forces during ∆T = 15 s is rather small and does only show up at lower speeds in the resulting distribution when longer sampling intervals are used.
In the second regime with high friction (η i = η med = 500 nN µm −2 s) (FIG. 6) the average speed with ∆T = 1 s and ∆T = 15 s match each other (not shown). The resulting speed distribution usually fit better to the experimental ones, in particular, because low speed contributions now are shifted to the peak resulting in a sharper speed distribution. Moreover the components of the speed distribution in the low speed range are smaller as the active forces can more easily overcome the JKR forces such that cells are less likely to be immobilized in an cellular environment of high density. However, the active forces are higher than measured experimentally [35,82,84]. Yet, the parameter set is computationally favorable as the high speed fluctuations are missing. This allows for larger time steps in the simulation.
The effect of active forces F act ij resulting from the constriction ring model is relevant to the speed distribution. Without these active inter-cellular forces the speed distribution shows only sharp peaks (FIG. 6 lower panel). Also when the speed distribution of each lymphocyte subset is fitted using a homogeneous cell aggregate the peaks in the speed distribution tend to be sharper than in experiment. This could be compensated with larger active cell interaction forces F act ij which however produces the wrong speed distribution when B and T cells interact which each other (not shown).

Stable follicle formation
The simple model of PLF formation is able to generate stable follicle sizes of a few 100 µm diameter with roughly 10 4 B cells. The follicle forms around the exit spots engulfing them almost completely (FIG. 7 (a)). The adjacent T zone is crescent-shaped but tends to form a closed shell around the B cells. The T zone is basically determined as non-B zone and T cells occupy the remaining space by random migration around the follicle. However, they do not diffuse freely in the whole space due to the chemoattraction to the exit spots. Thus this toy model of PLF formation generates a tissue in flow equilibrium with dynamically generated sources of cell attraction and keeps B and T cells separated.

Effect of receptor internalization
When studying a stable follicle including chemotaxis towards exit spots the effect of CXCR5 internalization dynamics on the follicle shape is not visible. In order to isolate the structural effect of CXCR5 receptor internalization, the chemoattraction to the exit spots is switched off. Then a flow equilibrium of B cells is no longer possible as the efflux of cells is strongly suppressed (not shown). This reflects that randomly migrating cells rarely find the exit spots. Consequently, cell influx and efflux are completely shut off which imposes a constant B cell number. This allows to study the FRC-FDC dynamics together with the internalization dynamics of CXCR5 alone. It is found that the internalization dynamics on its own is a source of instability for the system as the forming follicle cannot maintain its shape (FIG. 7 (b)).
Note, that this result is not affected by considering the internalization dynamics under B cell flow conditions. This can be investigated by compensating the low efflux of B cells, which arises due to switched off chemoattraction to the exit spots, by distribution of the exit spots over a large area. The instability is even stronger because the outflowing B cells act as an additional sink for CXCL13 by taking surface bound CXCL13 out of the tissue. The effect of receptor internalization can be understood by comparing the results to a simulation without CXCR5 internalization. If the binding state of the chemoattractant receptors is frozen the cells would migrate like being in a potential attracting them to FDC. This is rather independent of the B cell-induced dynamics between FDC and FRC because the migration is much faster such that the B cell distribution can equilibrate. Thus the distribution of B cells and FDC would match each other and the follicle would be stable. The internalization dynamics can however change the 'potential' formed by the chemoattractant on a time scale comparable to the cell migration (FIG. 7(c)). As the cells act as sinks for the chemoattractant they are always generating local gradients away from their current position. This cannot be overcome by diffusion of the chemoattractant which is slower than the cell migration on the corresponding cellular length scale. Note that the instability only occurs for a sufficiently large number of B cells because otherwise the sink for the chemoattractant is too small to change CXCL13 gradients significantly.
The instability of the follicle shape is amplified by the FRC-FDC dynamics. As the B cells spread out preferentially at the boundary of the follicle, where CXCL13 concentrations are rather low, they extent beyond the follicle border and increase the B cell density outside the follicle. Thus, the area covered by the B cells is bigger than expected from a dense packing of the B cells in the FDC area. This culminates in the generation of new FDC. As on one hand B cells are very motile and on the other hand FDC need some time to dedifferentiate to FRC, on long term this results in an FDC network which is bigger than the volume required by the number of B cells. Thus the whole PLF becomes unstable in shape following the changed 'center of mass' of the FDC network as determined by the concentration peaks of CXCL13.

Discussion
A model architecture was presented that allows to simulate fast cell migration taking into account subtle effects of chemoattractants and cell influx and efflux from tissues. Cells are represented as visco-elastic objects. They interact with adjacent cells by passive and active forces. The model allows to simulate detailed mechanics of single cells and individual coupling of internal cell dynamics to cell mechanics as well as to contact-and long-range-interactions. The use of an underlying regular triangulation permits a continuous description of cell positions and motility. It provides an efficient way of defining the neighborhood topology for cells of different size in tissue of different density. The advantage of the regular triangulation is to provide an effective method to implement fluxes of highly motile cells by dynamical modifications of the triangulation. In contrast, lattice-based architectures have difficulties to model fast cell migration. Either complicated rules are needed or a lattice-gas model is required [85]. In both cases the parameters coupled with the modifications are not easily converted to observable quantities. Using a physical representation of cells allows to directly access the involved parameters by experiment. This consequently requires a lattice-free description as implemented in the present model based on a regular triangulation. The triangulation is complemented by a grid for solubles like chemokines. The reaction-diffusion-equations are solved on that grid. Thus, the model acquires a hybrid model architecture. The model is based on parameters with physiological well-defined meaning such that the number of parameters is considerably reduced if the corresponding experimental values are known. This is, indeed, the case for most parameters in the present test-application.
The biomechanics of cellular aggregates of fast migrating cells has three components: passive mechanical interaction described by JKR forces, active forces exerted on extra-cellular matrix, and active forces directly acting between cells. In general the active forces dominate the JKR forces which allows cell movement in dense configurations. The speed distribution with high friction and consequently large active forces allows cells to more easily overcome the repulsion generated by JKR forces. Thus, the high friction regime mimics large deformations of cells without explicitely taking into account the shape of migrating lymphocytes which elongate in the direction of motion [26][27][28][29][30]. Alternatively, the elastic modulus E i of cells could be scaled to appropriate values below the physiologic values and correspondingly decreased adhesion parameters σ i . The constriction ring model takes into account lateral force generation by actively migrating cells which are ignored in other force-based models that consider only the generated net force [2,4]. It is shown that the cell interaction due to the exchange of active forces is the major determinant of the width of the speed distributions. Thus this interaction might be interpreted as a Brownian motion-like fluctuation of the cell speed also it is deterministic. The parameter coupled with the constriction ring is the pressure the pressure p * i (10). The force F act i (φ i ) exerted to the extra-cellular matrix may be either generated by the integrated pressure of the constriction ring model directly exerted to the matrix or by the filopodia dynamics of the classical three-step migration model (reviewed in [77,86]). Recent data demonstrates that both migration modes act in lymphocytes [87] although the contribution of the constriction ring is greater during chemotactic responses of migrating lymphocytes [87][88][89]. Especially, the characteristic time of the cytoskeleton dynamics of the ring coincides nicely with the persistence time [41,87]. The persistence time, the time interval during which the direction of active forces is constant, is likely to be an intrinsic feature of the underlying cytoskeleton dynamics [39,41,87,90,91].
The test-application for the full model involves a simple approach to the formation of PLF and has proven the functionality of the model. More specifically, the motility properties of the cells could be adjusted to in vivo motility data of lymphocytes in lymphoid follicles. The model generates homeostatic follicles in a flow equilibrium of reasonable size of few 100 µm [50,[92][93][94][95] demonstrating that the model architecture is able to describe systems that are in a cell flow equilibrium of fast migrating cells.
The effect of CXCR5 receptor internalization was demonstrated under non-flow conditions. Within the investigated parameter range (see appendix) the dynamics of CXCR5 and CXCL13 leads to a modification of the concentration profile such that cells temporarily migrate away from sources. This leads to a frequent change in shape of the cell aggregate. Such a behavior was not observed in vivo [28]. Therefore either internalization of CXCR5 does not occur in the PLF system or mechanisms exist that counter-balance the shape dynamics. Higher CXCL13 concentrations can overcome this effect as the internalization cannot act as a sufficient sink to locally invert the gradient of the chemoattractant (not shown). However, the cells then have strongly downregulated receptors such that they loose there chemotactic sensitivity, i.e. cells start to migrate randomly preventing efficient aggregation on long terms. We observed a stable follicle shape under flow conditions when the B cells chemotax to the exit spots without internalization dynamics. Thus, chemotaxis to the exit spots works as an independent attractor for B cells that can compenstate for locally inverted gradients of CXCL13 resulting from the receptor internalization of CXCR5.
The main purpose of this article was to establish a simulation platform suitable for modeling of homeostatic tissue dynamics on a cellular and subcellular level involving large numbers of fast migrating cells. The results have demonstrated that the proposed model design can cope with the complexity that occurs in tissue exhibiting a flow equilibrium of fast migrating cells. In particular, it clearly separates the various time and length scales and allows to localize the origin of emerging properties on the tissue level as well as on the cellular and molecular level. We, therefore, consider this simulation tool to be a suitable instrument for the analysis of morphogenesis of highly dynamic tissue. However, concerning the application to follicle formation, the assumed mechanisms are only subparts of full PLF formation dynamics that need to be refined in order to reproduce microanatomical data and to make realistic and quantitative predictions, which is left for future research.

Acknowledgments
FIAS is supported by the ALTANA AG. Meyer-Hermann is supported by the EU-NEST project MAMOCELL within FP6. *

A. Estimate parameters for chemoattractant reaction-diffusion system
The parameters used in the ODE system (15) are not known but can be estimated from data of similar systems. The dissociation constant K d for chemoattractants and their receptors are mostly measured for other chemoattractants than the ones used here (CCL21 and CXCL13) ( Table 2). The values for K d range from 0.2 nM to 5 nM [96][97][98][99]. The dissociation constant for CCR7 with CCL21 has been measured to be 1.6 nM [100]. Considering the range of these values it is likely that K d adopts values in this range in the present system. A less favorable situation exists for the reaction rate constant k on (k off = k on K d ). Only few data exists which spread over several orders of magnitude. Values for the association rate k on are available for CXCL12 binding to fibronectin [98] (2.5 · 10 5 M −1 s −1 ) . The off rate k off measured for CXCL12 binding to fibronectin is 6.5 · 10 −3 s −1 . The two rate constants k i and k r associated with receptor internalization and recycling (see (15)) can be estimated from experimental data on receptor desensitization and resensitization experiments which are reviewed recently [45]. It is assumed that all non-internalized receptors have bound the ligand when high chemoattractant concentrations far above K d are used as done in the experiment. The only equation that is left for the desensitization process is then (15) (R = 0) with the solution The internalized receptor fraction r * at equilibrium becomes   (15).
Assuming that the internalization process is not in steady state -and solving (17) numericallydoesn't change the results very much compared to the experimental uncertainty. The numerical results of PLF-formation are not sensitive to these parameters (data not shown). R tot is not known explicitly. From similar receptors the number of CCR7 molecules on T cells has been estimated to be 10 5 per cell [100] and 10 4 for B cells as indicated by the studies that find a factor 10 difference of CCR7 levels between B and T cells [101]. Note that the values cited above have to be converted into the receptor concentration R tot based on the cell densities of the corresponding lymphocyte type.