Adhesion and volume constraints via nonlocal interactions determine cell organisation and migration proﬁles

The description of the cell spatial pattern and characteristic distances is fundamental in a wide range of physio-pathological biological phenomena, from morphogenesis to cancer growth. Discrete particle models are widely used in this ﬁeld, since they are focused on the cell-level of abstraction and are able to preserve the identity of single individuals reproducing their behavior. In particular, a fundamental role in determining the usefulness and the realism of a particle mathematical approach is played by the choice of the intercellular pairwise interaction kernel and by the estimate of its parameters. The aim of the paper is to demonstrate how the concept of H-stability, deriving from statistical mechanics, can have im- portant implications in this respect. For any given interaction kernel, it in fact allows to a priori predict the regions of the free parameter space that result in stable conﬁgurations of the system characterized by a ﬁnite and strictly positive minimal interparticle distance, which is fundamental when dealing with biological phenomena. The proposed analytical arguments are indeed able to restrict the range of possible variations of selected model coeﬃcients, whose exact estimate however requires further investigations (e.g., ﬁtting with empirical data), as illustrated in this paper by series of representative simulations deal- ing with cell colony reorganization, sorting phenomena and zebraﬁsh embryonic development.


Introduction
An accurate description of the spatial pattern of cell aggregates is a fundamental issue in theoretical biology. The spatial cell configuration and characteristic distances are in fact at the basis of a wide range of biological processes, i.e., from morphogenesis to cancer growth and invasion. For example, defects in the spatial organization of multipotent stem cells in animal embryos lead to severe malformations of adult organs ( Ilina and Friedl, 2009 ). Further, the compact configuration of epithelial monolayers is fundamental in wound healing scenarios. Finally, the dispersion of highly motile malignant individuals triggers the metastatic transition of tumor progression .
From a mathematical point of view, the spatial distribution of cells forming aggregates can be well approximated by discrete models, which actually approach the biological problem focusing on the cell-level of abstraction and preserve the identity and the behavior of individual elements (for comprehensive reviews the * Corresponding author. Tel.: +393471827593. E-mail addresses: carrillo@imperial.ac.uk (J.A. Carrillo), annachiara.colombi @polito.it (A. Colombi), marco.scianna@polito.it , marcosci1@alice.it (M. Scianna). reader is referred to Alber et al. (2003) ; Anderson and Rejniak (2007) ; Deutsch et al. (2007) ; Drasdo (2003b) ). In more details, these techniques represent biological elements as one or a set of discrete units, being individual morphology restricted according to some underlying assumptions. Among discrete approaches, we here focus on particle-based models, where the biological individual is represented by a material point with concentrated mass and identified by its position in space, in contrast to methods that allow for a description of the cell membrane and its morphology as vertex-based models ( Fletcher et al., 2014 ) or Cellular Potts Models (CPM, see Balter et al. (2007) ; Glazier et al. (2007) ). For particle or agent-based models, the cells move according to ordinary differential equations (ODEs), with the phenomenological postulation of either acceleration (second-order models) or velocity (first-order models) contributions. In both cases, among the possible migratory components that can be included, one of the main important model ingredients is the term relative to direct intercellular pairwise interactions that can be described by a proper kernel (potential). This contribution in cell behavior is typically the combination of adhesive/repulsive mechanisms and can lead to unrealistic cell collapse or dispersion or more realistic patterns with optimal individual spacing.
The definition of proper interaction kernels, as well as an accurate estimate of the relative parameters, is therefore mandatory when developing a particle model for biological problems. In particular, the issue relative to the parameter estimate is very relevant in theoretical biology. In fact, regardless of the specific implemented approach, a direct one-to-one correspondence between all the model parameters and experimental quantities is not straightforward (as commented also in Balter et al. (2007) ; Glazier et al. (2007) ; Merks and Koolwijk (2009) , in the case of CPM). The different model coefficients often interfere with each other in an intricate way, and therefore simultaneous parameter fittings are typically needed, which can be done with large-scale massive preliminary simulations.
Objective of the work. The main objective of this paper is indeed to propose a procedure that improves the strategy to choose proper adhesive/repulsive kernels by identifying a physically admissible subset of the free parameter space which gives rise to realistic system configurations in terms of inter-agent spacing. In order to do this, we will take advantage of the concept of H-stability of particle interaction kernels and potentials, defined in the statistical mechanics ( Ruelle, 1969 ), and already used in the modeling of swarming and collective migration of animal population ( Carrillo et al., 2010b;Kolokolnikov et al., 2013 ). The determination of the H-stability properties allows in fact to predict the stable configurations, in terms of particle spacing, of the system of interest. Of particular relevance is the fact that the H-stability condition of interaction kernels translates into a constraint involving the values of the relative interaction parameters, therefore reducing the range of their possible variations if we want to obtain realistic patterns. It is useful to clarify that the proposed analytical results do not shed light either on the dynamics of the aggregate or on the exact position of the single agents at the equilibrium.
Structure of the work. The rest of the paper is then organized as it follows. In Section 2 , we will present a general first-order particle model (formally derived by a second-order approach) and discuss some possible velocity contributions that can be included in the modeling framework. In Section 3 , we will focus on the interaction velocity term. In this respect, we will propose a class of interaction kernels/potentials, commenting the underlying biological hypothesis and, in Section 3.2 , we will introduce the concept of H-stability of the system and its implications in the classification (in terms of agent spacing) of the stable configurations of a particle system and on the estimate of the interaction parameters. In Section 3.3 , such analytical results will be supported by proper sets of simulations. Section 3.4 will instead show that, once restricted the free space of the interaction model parameters, a more detailed coefficient estimate can only result by a data fitting with experimental/biological quantities. Finally, in order to not remain on a pure conceptual level, we will finally show in Section 4 that the proposed analytical procedure can be important also in the case of more realistic biological applications: for instance, we will reproduce the early migration of the zebrafish lateral line primordium, whose dynamics are fundamental for the correct embryonic development of the animal, and show how the use of H-stable interaction potentials is crucial to reproduce typical migration patterns.

Basic Mathematical Model
We start by considering a biological system composed by N cells of the same phenotype/lineage, i.e., characterized by the same biophysical properties (mass and dimension) and behavior. We anticipate that in the following sections, we will extend such a model framework in the case of multiple differentiated cell populations. As previously introduced, each individual is here represented as a discrete entity, i.e., as a material particle with concentrated mass, say m , and characterized by its position in space, x i (t) ∈ R 2 with i = 1 , . . . , N, assuming a planar cell distribution. The configuration of the overall system at a given instant t can therefore be given by the vector: To approach the dynamics of a generic cell i , we start from a general second-order particle model: (1) The second and the third terms on the left hand side of Eq. (1) describe damping mechanisms related to the friction forces resulting from the contact of cell i with, respectively, the substrate and the other individuals, as λ cs i and λ cc i j are the corresponding coefficients (which may be, for instance, constant parameters in common for all cells ( Costanzo et al., 2015 ) or time-dependent individuallyspecific tensors ( Liedekerke et al., 2015 )). These contributions have the effect to slow down individual movement and eventually to increase the characteristic time scale of the overall cell collective patterning. On the right hand side, F i instead denotes the sum of all forces influencing cell behavior. However, in order to simplify the picture, we can first notice that cells move in extremely viscous environments, characterized by very small Reynolds numbers: inertial effects in cell dynamics can therefore be neglected, if a sufficiently large observation time is considered ( Liedekerke et al., 2015;Odell et al., 1981 ). In fact, in these conditions, biological cells can maintain a persistent ballistic locomotion only for a substantially small time, giving rise to straight displacements shorter than their typical dimensions. These considerations allow to drop the inertial term in (1) and to employ a first-order model, where the velocity of an individual, and not its acceleration, is proportional to the acting forces. Such a relation, called overdamped force-velocity response , is at the basis of a number of other discrete/IBM approaches (see Drasdo (2003b) ; Scianna and Preziosi (2012) and references therein for comments). As seen, cell-cell friction is in principle an important contribution: however, such a damping term is of particular relevance in the case of threedimensional settings, where the adhesive surface between cells is significantly large ( Drasdo, 2003b ). In fact, in planar domains, as those considered in this work, cell membranes are instead mainly in contact with extracellular elements, e.g., flat matrix substrates. It is further not too restrictive to assume, at least in a first approximation and in a more conceptual work, that the relative velocity between pairs of cells is negligible with respect to cell-substrate friction. For these reasons, we opt to neglect also the third term on the left hand side of Eq. (1) , thereby obtaining the following law of cell motion: ( Eq. (2) implies indeed that cell dynamics can be described by a direct phenomenological postulation of the velocity contributions, which have to take into account of the cell-substrate friction coefficient, possibly included within their characteristic parameters. It is finally useful to remark that, once cells move with constant velocity, the shape of their migration profile is given by the stationary states of the model (2) : they derive from the balance of forces acting on each individual (i.e., F i (t) = 0 for all i ) and are exactly the same that would result by the corresponding second-order approach. We now assume that the velocity of each cell i results from the superimposition of different, biologically reasonable, contributions: (3) In Eq. (3) , v dir i is a directional component in individual motion, that may reflect chemotactic, haptotactic, or durotactic mechanisms (i.e., cell locomotion towards substrate regions with increasing concentrations/density of soluble or insoluble chemicals, or of rigid extracellular matrix components, respectively). v int i models instead individual dynamics resulting from direct intercellular interactions, which are the main topic of this article and will be dealt in more details below. v rand i is a noise term, that takes into account that biological elements (not only cells but also bacteria and other organisms) explore their surrounding environment and crawl in a random fashion. Finally, v other i possibly includes other cell migratory mechanisms, such as Cucker-Smale type alignment contributions ( Cucker and Smale, 2007 ), employed for instance to describe dynamics of swarms ( Carrillo et al., 2010a ), epithelial fronts ( Sepúlveda et al., 2013 ) and zebrafish embryos ( Costanzo et al., 2015 ), although in second-order approaches.

Intercellular Interaction Velocity Contribution
Of the above-introduced spectrum of cell velocity components, let us now focus on the term relative to cell-cell direct interaction instances, i.e., v int . A detailed analysis of its structure and effects, as well as a reasonable estimate of the relative characteristic parameters, is in fact critical in the mathematical description of most biological phenomena, since cell-cell interactions are at the basis of system configuration, patterning, and development ( Friedl and Bröcker, 20 0 0 ).
More precisely, we assume that cell behavior due to direct intercellular interactions results from the superposition of the contributions of pairwise forces, which consistently depend on the relative distance between the two individuals involved and result aligned to line ideally connecting them. Based on these working hypothesis, the cell-cell interaction velocity term can be given by: where | · | identifies the Euclidean norm and K : R + → R (that has units μm/s) is an interaction kernel. In (4) , we consistently avoid cell self-interactions as well.
First-order models with nonlocal attractive and repulsive terms to include cell adhesion and volume effects have been already proposed in the literature of cell interactions for cancer invasion models ( Colombi et al., 2015a( Colombi et al., , 2017Domschke et al., 2014;Gerisch and Chaplain, 2008;Painter et al., 2015 ), zebrafish lateral line patterning ( Volkening and Sandstede, 2015 ), and cell sorting in heterogenous cell populations ( Burger et al., 2017;Carrillo et al., 2017;Murakawa and Togashi, 2015 ). Some of this works do not deal with agent-based models but with their continuum macroscopic limit or even with an hybridization of both description levels ( Colombi et al., 2015a( Colombi et al., , 2017. For instance, at the macroscopic level, repulsion is taken into account by (nonlinear) diffusion or drift saturation terms ( Burger et al., 2017;Calvez and Carrillo, 2006;Carrillo et al., 2017;Hillen and Painter, 2001;Painter and Hillen, 2002;Shigesada et al., 1979 ), which can be obtained from particle-based nonlocal repulsive models in the right scaling limit ( Bodnar and Velázquez, 2013;Oelschläger, 1990 ).

Definition of the Intercellular Interaction Kernel
Despite cell-cell direct forces may result from several mechanisms, we hereafter take into account only repulsive interactions , which reproduce cell resistance to compression due to size and mechanical arguments Serini et al., 2003 ), and attractive interactions , which conversely implement cell-cell adhesiveness, that relies upon the expression and the activity of transmembrane molecules (i.e., cadherins). In this respect, denoting as x, y the position in space of a pair of particles, we say that K implements a cell repulsive behavior when K(| x − y | ) < 0 and an attractive one when K(| x − y | ) > 0 , see Fig. 1 (right panel). Coherently, cell resistance to compression enters the picture if the relative distance between the interacting individuals is lower than a minimal vital space, i.e., lower than a mean cell diameter, hereafter denoted by d r (see Fig. 1 , left panels). On the other hand, cell-cell adhesiveness is active if the relative distance between two interacting particles is large enough to avoid cell-cell repulsion, i.e., | x − y | > d r , but sufficiently small to allow the formation of bonds between cell membrane protrusions, i.e., smaller than the maximal possible extension of cell deformable adhesive structures, which is hereafter denoted by d a (see Fig. 1 , central panels). Finally, when two individuals fall too apart one to each other, i.e., their distance is larger than d a , they do not interact: in this case, K(| x − y | ) is set equal to 0, see Fig. 1 , right panel.
In principle, there are many possible choices for the explicit form of the interaction kernel K introduced in Eq. (4) : however, they have to reasonably and accurately describe the biological phenomenon of interest. In particular, we hereafter consider a set of interaction kernels characterized by a parabolic shape in the positive/attractive part and a differentiated behavior in the negative/repulsive part ( Colombi et al., 2017( Colombi et al., , 2015b. Specifically, the cell repulsive neighborhood accounts for distinct compressibility of cell nucleus and cytoplasm, with the former that is typically less squeezable than the latter, that is conversely highly deformable, as widely proven in the experimental literature ( Ilina and Friedl, 2009;Wolf et al., 2003 ). In this respect, we denote the nucleus diameter by d n (that is reasonably close to cell radius, i.e., d n = d r / 2 ) and assume a different shape of the repulsive part of interaction kernel depending on whether the intercellular distance is lower than d n or falls in the range [ d n , d r ].
According to all the previous considerations, we indeed hereafter deal with the following family of interaction kernels: The functions defined in Eq. (5) , and plotted in Fig. 2 (left panel), are intrinsically multiparametric, since they are characterized by the following set of coefficients: To decrease the complexity of the problem, we can first restrict the dimension of the free parameter space with phenomenological arguments and observations. The values relative to cell dimensions, i.e., d r and d a , can be in fact easily taken from the experimental literature according to the phenomenon of interest and to the specific type of cells involved. This reduces the set of free parameters to In this respect, F r and F a (that have both units μm/s) are related to the repulsion and adhesive strengths of cell-cell interactions, Left panels: Repulsive interactions affect the dynamics of two cells when their relative distance is lower than a mean cell body diameter d r . The repulsive neighborhood of each cell i , with i = 1 , . . . , N, then consists of a round area with radius d r centered at the actual position of particle i . Central panels: Attractive interactions arise between cells whose relative distance is larger than the mean cell body diameter d r and lower than the maximal possible extension of cell membrane adhesive structures d a .
The cell adhesive neighborhood then consists of a circular ring with inner radius d r and external radius d a . Right panel: The radial interaction kernel K is assumed negative when repulsive, i.e., K(| x − y | ) < 0 when | x − y | < d r , and positive when attractive, i.e., Finally, cells whose relative distance is too large do not interact: consistently, we set

Fig. 2.
Representation of the family of intercellular pairwise interaction kernels K ( r ) defined in Eq. (5) and the corresponding potentials u ( r ), see Appendix A . All the proposed functions present a parabolic behavior in the attractive part and a differentiated shape in the repulsive part. In particular, the repulsive trend near the origin (i.e., when the distance between the two interacting cells is lower than d n = d r / 2 ) is established by the parameter s ∈ (0, 2]. In the figure, we show the shape of some representative kernels/potentials obtained by setting s = 2 , 1 . 75 , 1 . 5 , 1 . 25 , 1 , 0 . 75 , 0 . 5 , 0 . 25 , we focus on in the following sections. respectively, being also scaled by the cell-substrate friction coefficient (see before). In more details, F r is proportional to an intrinsic cell incompressibility/stiffness, that is necessary to allow the single individuals to realistically preserve their dimensions and to avoid overlaps; F a is instead correlated to the cell-cell adhesiveness which, from a molecular viewpoint, is locally regulated by the expression and the activation of specific cadherin molecules. Finally, the parameter s characterizes the kernel behavior near the origin (see Fig. 2 ), being broadly related to the nucleus stiffness. The values of F r , F a and s do not indeed have a clear and direct biological counter part and are therefore difficult to estimate, although they have to result in realistic cell velocities. However, an analysis based on arguments deriving from statistical mechanics is able to a priori predict the regions of remaining free parameter space that give rise to stable system configurations characterized by optimal cell spacing, fundamental when modeling of biological phenomena. The proposed study indeed facilitate the parameter esti-mate, at least by reducing the possible variations of selected model coefficients.

H-Stabilty of the Intercellular Interaction Kernel: Analytical Results
With this aim, and in order to simplify, let us start by neglecting all the other velocity contributions, except from the one deriving from direct intercellular interactions. We remark that the inclusion of other cell behavior is possible, as it will be shown in the following sections. Taking into account of (4) , Eq. (2) can be rewritten in the following form To discuss some the theoretical results concerning the stationary states of problem (6) , notice that the family of interaction kernels K introduced in Eq. (5) can always be considered as being derived from a scalar interaction potential u : R + → R such that u (r) = K(r) and that it is possible to define a vector potential U : for all i = 1 , ...N. The system of equations (7) has the structure of a gradient flow of the total potential energy In particular, E N is a Liapunov functional for (7) and therefore the stable stationary states of the problem are among (local or global) minimizers of the interaction energy. Let us now discuss some of the qualitative properties known for these (local) minimizers.
The first important theoretical result is that the behavior of the solution at the equilibrium is regulated by the repulsive part of the interactions, i.e., by the singularity of the interaction potential/kernel at the origin. This fact was studied in Balagué et al. (2013) ; Bertozzi et al. (2015) where it was shown that, as the potential gets more and more repulsive at the origin, the particles distribute in larger and larger regions. In other words, while mild repulsion may allow for clustering of particles, singular repulsion leads to regular distributions of particles in the plane and a well defined minimum interparticle distance.
The second theoretical aspect is relative to the concept of Hstable potentials, introduced in statistical mechanics, see Ruelle (1969) , and intimately related to the emergence of crystal-like behavior in ensembles of interaction particles. Assume u is a potential essentially negligible for large distances, i.e., such that lim r→ + ∞ u (r) = 0 : if u is H-stable, we have that at the equilibrium the minimal interparticle distance is bounded below by a finite fixed positive value regardless of the total amount of individuals N . On the opposite, if the potential u is not H-stable, sometimes also called catastrophic in the statistical mechanics literature, the minimal interparticle distance at the equilibrium collapses to 0 as N → ∞ . An easy-to-check criterium to detect the H-stability of a potential is given in Ruelle (1969) and summarized in the next Theorem, see also Carrillo et al. (2016) ;nizo et al. (2015) ; nizo and Patacchini (2016) for further discussions and results in the not Hstable cases.
Basing on this characterization, we can state the following:

Corollary 2. The family of interaction potentials defined in Eq. (5) is
For the sake of completeness the proof of the Corollary is given in Appendix A .
Taking into account that, as seen, the values of the interaction radii d r and d a are defined according to cell phenotype, the above analytical result states that the H-stability (and therefore the overall system behavior) translates into a constraint on the ratio between the repulsive and adhesive interactions strengths (i.e., F r and F a ), fixed the value of s . In other words, for any given interaction kernel, only pairs of values F r and F a satisfying (8) will give a stable system configuration characterized by a finite and strictly positive interparticle distance, i.e., a biologically realistic model.

H-Stability of the Intercellular Interaction Kernel: Computational Results
Let us devote the rest of this section to computationally support our previous analytical considerations, i.e., to showcase how the final pattern of pairwise interacting particles is governed by the H-stability of the interaction kernel and its behavior at the origin. To do this, we hereafter perform numerical tests that reproduce the evolution of aggregates of N component cells, whose dynamics are regulated by model (7) with interaction kernels K defined as in Eq. (5) . In more details, our strategy is to vary the values of the free interaction parameters, i.e., s, F r , and F a and to explore the discrete stationary states of N -cell populations. In particular, the system behavior will be classified according to the following quantities which will be evaluated at the beginning of the simulation, i.e., at t = 0 (initial condition) and after t F , i.e., an observation time sufficiently large to allow the component cells to reach a stable equilibrium configuration. We remark that the measures introduced in (9) -(10) have a well-defined biological meaning: d min is the minimal intercellular distance, whereas d max gives in fact the extension of the overall particle distribution. If its value is above a threshold value (consistently close to the nucleus diameter d n ), all cells have a sufficient space to remain viable and survive. Otherwise, they can be considered overlapped or dramatically compressed and therefore suppressed.
Variations of N in the cases of different values of s and F r / F a . We first investigate the effect on the equilibrium configuration of cell colonies of variations in the overall number of agents in the case of either H-stable or not interaction kernels K . A series of numerical realizations is indeed performed on systems formed by N = 50 , 100 , 200 particles which evolves following (6) with interaction kernel defined in (5) . Regardless of their number, the agents are initially distributed in an almost round area of radius equal to 60 μm (i.e., d max (0) = 120 μm), see Fig. 3 . In particular, in all cases, the initial configuration of the aggregates is such that for each cell i = 1 , . . . , N the minimal interparticle distance d min (0), defined in Eq. (9) , is not zero (in order to avoid that distinct individuals are initially located at the same position) and lower than d a (so that each cell is initially able to interact at least with another cell). According to Colombi et al. (2015aColombi et al. ( , 2017 and biological references therein, the cell typical dimensions are finally taken as d r = 20 μm (so that d n = 10 μm ) and d a = 60 μm .
Specifically, we focus only on three possible choices for the behavior of the interaction kernel near the origin, determined by s = 1 . 75 (more regular than hyperbolic, which results in F * = 37 , 26 ), s = 1 (hyperbolic, which results in F * = 31 , 32 ), and s = 0 . 25 (more area of radius equal to 60 μm (i.e., d max (0) = 120 μm ). The initial minimal intercellular distance is not null to avoid overlapping and such that each cell is able to interact at least with another individual, i.e., d max (0) ∈ (0, d a ].

Fig. 4.
Equilibrium configuration of cell aggregates formed by an increasing number of component cells, observed for selected representative choices of the shape of the interaction kernel at the origin (i.e., s = 1 . 75 , 1 , 0 . 25 ) and of the ratio between the interaction strengths parameters (i.e., F r / F a , which makes the interaction kernels H-stable or not). In the sake of completeness, in each panel, we indicate the stable values of d max ( t F ) and d min ( t F ). singular than hyperbolic, which results in F * = 15 , 36 ). Simultaneously, we consider two representative values of the ratio F r / F a , i.e., F r /F a = 10 and 100, that respectively makes the considered interaction kernels not H-stable and H-stable.
The computational results, summarized in Fig. 4 , show that in the case of not H-stable kernels: -the minimal intracellular distance at the equilibrium, i.e., d min ( t F ), dramatically decreases upon increments in N ; -the overall extension of the cell aggregate, i.e., d max ( t F ), remains almost constant regardless of the number of particles forming the population.
In more details, for any fixed value of N , if the interaction kernel is not H-stable, the individuals tend to collapse, towards a colony of d max ( t F ) ≈ 70 μm. Slightly differences can be however observed due to the characterization of K near the origin. In particular, we have that, for any tested N , by setting s = 1 . 75 , cells organize into several compressed clusters, with consistent cell overlapping, i.e., d min ( t F ) < 2 μm. A more homogeneous spatial distribution of cells at the equilibrium is then observed for lower values of s (i.e., for interaction kernels more singular near the origin): however, also in these cases, the stable intercellular distance is still too small to have a realistic individual survival ( d min ( t F ) < d n ). From an applicative/biological perspective, this system behavior corresponds, for instance, to the fact that if we culture more and more cells within the same, say, Petri dish, they would organize into increasingly denser colonies, collapsing and overlapping one from each other. That would be obviously an unrealistical situation.
On the opposite, in the case of H-stable particle systems, increments in N result in increments in the overall dimension of the cell aggregate at the equilibrium (i.e., of d max ( t F )) whereas the characteristic intercellular minimal distance d min ( t F ) does not vary: it remains close to the repulsive radius d r , being independent from N . In particular, in all H-stable cases, we have cell colony re-organization and enlargement (i.e., d max ( t F ) > d max (0)), with the stable configuration characterized by a minimum intercellular dis-  5) ). In particular, the adhesive coefficient F a is constantly assumed equal to 1 μm/s, so that the value of the ratio F r / F a results equal to the value of the repulsive parameter F r . According to Corollary 2 , for each value of s , the interaction kernel/potential is not H-stable when F r / F a < F * and H-stable when F r / F a > F * . In this respect, the triangles indicate for each curve (i.e., for each s ) the projection of corresponding value of F * . We remark that values of d min ( t F ) lower than d n results in unrealistic cell overlapping.
tance that allows individuals to have a sufficient vital space, which does not significantly change upon variations of s . Biologically, it is a reasonable behavior, since, if we add more cells, then an aggregate tends to enlarge and to occupy extended parts of the substrate to allow each individual to have a sufficient vital space.
Variations of s and F r / F a . Now, we consider the representative cell aggregate constituted by N = 100 individuals, with the initial configuration defined in the previous section and shown in Fig. 4 , middle panel, and investigate its stable configuration for a wider range of ratios F r / F a , in the case of the following values of s (and of F * which, for any given s , defines the minimum ratio between the repulsive and adhesive interactions strengths that makes the interaction kernel K H-stable, see Eq. (8)  To avoid overcomplications, in all realizations, we keep the strength of cell-cell adhesiveness F a constantly equal to 1 μm/s and vary the value of the intrinsic cell resistance to compression F r . This means that we are fixing the adhesive characteristics of cell-cell interactions, focusing on the repulsive part. This is not a shortcoming, since the discriminating quantity is the ratio F r / F a . In particular, we test the following values: F r = 1 , 10 , 20 , 30 , 40 , 50 , 100 , 10 0 0 μm / s .
In agreement with the above theoretical and computational considerations, the graphs in Fig. 5 clearly show that the Hstability of the interaction kernel K determines the characteristic measures of the final configuration of the system. In fact, for any fixed value of s , a biologically realistic crystalline-like cell pattern, where the component individuals enlarge and stabilize at a sufficiently large distance without collapsing, can be obtained only for F r / F a > F * , i.e., if K is H-stable. In particular, in all these cases, d min (0) < d n < d min ( t F ) < d r , i.e., the minimal intercellular distance at equilibrium is larger than cell nuclear dimensions and, when F r / F a is sufficiently higher than F * , close to the repulsion radius d r (the cell body diameter), regardless of the value of s . From a biological view point, this means that the aggregate stabilizes in a colony where the component cells have a sufficient vital space. On the opposite, for all values of s , ratios of F r / F a lower than F * lead to cell aggregate dramatic collapse, as we constantly have d min ( t F ) < d min (0) < d n , with d min ( t F ) close to 0 for F r / F a substantially low. The above described system behavior is confirmed by the value of d max ( t F ) as well. The overall aggregate diameter is in fact much larger in the cases of H-stable interaction kernels. It is finally useful to remark that the simulation results are not dependent on the exact initial distribution of cells, satisfying the conditions on the initial interparticle distance (i.e., d min (0) ∈ (0, d a ]).
The different sets of numerical realizations presented in this section allow to claim that for any given interaction kernel, only sets of parameters satisfying the H-stability criterium results in realistic (in terms of particle spacing) stable cell configurations. This is of some help in reducing the possible variations of the model interaction coefficients.

H-stable systems and Parameter Estimate
Once restricted with analytical arguments the space of the possible variations of the free coefficients relative to adhesive/repulsive cell behavior, a more detailed parameter estimate can be obtained only by focusing on a specific application and comparing critical quantities characterizing the virtual cell system with the proper experimental counterparts.
In this respect, let us deal with a quiescent spheroid of 817 ovarian cancer cells (OVCAR:39) plated on a two dimensional Petri dish, see Fig. 6 , where the cell nuclei have been colored in yellow. The aim of the following study it to find a parameter setting that allows a computational aggregate formed by the same number of individuals, which evolves following (7) with interaction kernel K belonging to the family defined in Eq. (5) , to stabilize in a configuration reasonably close to the experimental one. As done in the previous section, we take again d r = 20 μm , d a = 60 μm (these measures are also consistent with ovarian malignant cells, see Shield et al. (2009) ) and F a = 1 μm / s , so that the free coefficients again remain the pair ( F r , s ). The computational particle are initially randomly disposed in a round area of radius 120 μm, with d min (0) ∈ (0, d a ] (see the comments in Section 3.3 ).
We then compare the stable particle configurations resulting for different parameter values, which however satisfy the H-stability criterium of Corollary 2 , with the experimental colony by using the following critical quantities: where d min ( t F ) and d max ( t F ) are defined in Eqs. (9) -(10) . d ex min and d ex max denote instead the minimal and maximal intercellular distance measured in the biological system (by a post-processing technique able to map position and center of mass of the cell nuclei, see Fig. 6  > 0, then the numerical results overestimate the real minimal intercellular distances (rsp., overall aggregate extension). On the opposite, if d min (rsp., d max ) < 0, the corresponding experimental value is underestimated. Moreover, if | d min | and | d max | are both smaller than 1, the discrepancy between the numerical configuration and experimental one is lower than a cell diameter. The values of d min and d max obtained from selected parameter sets are finally summarized in Fig. 7 . It is first possible to observe that, in all cases, the experimental minimal intercellular distance is slightly overestimated, i.e., 0.3 < d min < 0.6, for any choice of s and F r . In particular, in accordance with the results proposed in Section 3.2 , d min is not dependent on the behavior of the kernel near the origin (i.e., on s , which as seen is a sort of a measure of the nucleus compressibility) and slowly increases upon increments in F r / F a .
On the opposite, the discrepancy of the overall extension between the experimental and the virtual colony significantly varies according to the specific set of parameters. In particular, regardless the value of s , we have an overestimation of the aggregate diameter for F r / F a > 200, with increment in the error upon increments in F r / F a . An increasing underestimation of the colony diameter is instead observed upon decrements in the ratio F r / F a starting from the value of F r /F a = 75 . Interestingly, in the case of F r /F a = 100 we have both types of error depending on the value of s , although their are acceptable (i.e., | d max | < 0.75 for any tested s ).
Taking all these results together, we can state that there is not a set of parameters that simultaneously minimizes both d min and d max . The minimal intercellular distance of the experimental population is in fact best approximated with s = 0 . 5 and F r /F a = 50 , while the overall extension of the OVCAR colony by the pair of values s = 1 and F r /F a = 100 .
The proposed series of simulations shows that the H-stability theory is necessary to predict the realism of stable cell configu-rations, but not sufficient to derive a specific parameter setting. In particular, in this case, even a data fitting with empirical quantities is not completely satisfactory since distinct groups of interaction coefficients optimize distinct experimental measures. In similar situations, the choice of the values to be given to the free model coefficients can be however suggested by the problem of interest. If, as in this case, we are dealing with tumor growth, the most critical measure (i.e., the one that has to be approximated in more details) is typically the overall extension of the lesion: this implies to opt for the parameter setting that gives the minimal error d max . Otherwise, if we are deling, for instance, with epithelial sheets for biomedical scaffolds in the case of wound healing, we may be instead interested in the compactness of cell aggregates, i.e., in the approximation of the minimal intercellular distance and therefore in minimizing d min .

Biological applications
This section will be devoted to show how the application of our theoretical procedure can help in modeling more complex biological phenomena, if it is supported by proper empirical considerations and data comparison. One one hand, we will deal with cell sorting phenomena, where different cell types are involved. On the other hand, we will analyse the embryonic development of the zebrafish posterior lateral line primordium, where the main goal will be to capture the migration of two types of mesenchymal cells, i.e., leader and follower individuals, that can crawl with different speeds along the myoseptum.

Cell sorting
In multicellular organisms, the relative adhesion of various cell types one to another or to noncellular components surrounding them is also fundamental. From the late 1950s, it has been widely noticed that during embryonic development the behavior of cell aggregates resembles that of viscous fluid. A random mixture of two types of embryonic cells, in fact, spontaneously reorganizes to reestablish coherent homogeneous tissues ( Armstrong and Parenti, 1972;Steinberg and Takeichi, 1994 ). A similar process is a key step also in the regeneration of normal animal from aggregates of dissociated cells of adult hydra ( Gierer et al., 1972 ). It also explains the layered structure of the embryonic retina. These phenomena,   8. Initial distribution of the two-population aggregate (i.e., composed of light "L" and dark "D" particles) employed to reproduce cell sorting phenomena. The aggregate is constituted by 800 cells divided in two groups equal in number, represented by circles of radius d r . The initial condition has been randomly generated such that the minimal intercellular distance at the initial instant is greater than 0 (i.e., cell overlapping is initially avoided) and lower than the maximal extension of cell filopodia. commonly called cell sorting , involve neither cell division nor differentiation, but are entirely caused by spatial rearrangements of cell positions due to differences in the specific adhesiveness, see also Graner and Glazier (1992) ; Steinberg and Takeichi (1994) and references therein. In particular, specific hierarchies of adhesive strengths lead to specific configurations of the cellular aggregate.
A simple and intuitive simulation reproducing biological cell sorting deals with a particle aggregate formed by two types of individuals, namely, light "L" and dark "D" (that hereafter will be graphically represented by white and black circles, respectively, see Fig. 8 ). The dynamics of each individual (regardless of its type) are then assumed to be entirely due to adhesive/repulsive interactions. In this respect, the evolution in time of the spatial distribution of the particle system is given by the extension to two populations of model (4) : . . . , N L and h = 1 , . . . , N D denote the actual positions of the "L" and "D" cells, respectively; while the interaction kernels K pq : R + → R , with p, q ∈ {L,D}, define how a cell of phenotype p reacts to the presence of a cell of phenotype q . The absence of other mechanisms involved in cell behavior is consistent with previous works on cell sorting based on different types of mathematical approaches ( Graner and Glazier, 1992;Mackey et al., 2014;Murakawa and Togashi, 2015 ). Each kernel K pq introduced in (12) where the coefficients F pq a are a measure of homotypic, i.e., for p = q, or heterotypic, i.e., for p = q , cell-cell adhesiveness. Fixed cell dimensions according to the problem of interest, the possible variations of the remaining free parameters, s, F r , F LL a , F DL a = F LD a (due to reciprocity of interactions) and F DD a can be further reduced following the analytical and numerical analysis proposed in the previous sections. To have biologically realistic stable configurations, the cell system has to satisfy the H-stability condition. This translates into a series of constraints between the interaction parameters: for each pair ( pq ) ∈ {LL, DL, DD} we have in fact that the ratio F r /F pq a has to be lower than the threshold F * deriving from the choice of the value of s (cf. Corollary 2 ).
Taking all these considerations into account, it is now possible to show how the H-stability theory is crucial to assure realistic cell spacing, but the exact individual distribution and pattern is completely determined by the specific parameter setting. In this respect, we start with a population formed by 800 cells, equally distributed between dark and light individuals (i.e., N L = N D = 400 ). They are initially placed within a round area with initial interparticle distances that avoid overlaps and allow each individual to interact at least with another one. Given again d a = 60 μm and d r = 20 μm, this time we fix the repulsive coefficient F r = 1 μm/s and vary The numerical results collected in the left column of Fig. 9 show that a proper interagent spacing is consistent in all realizations. Moreover, if the heterotypic adhesiveness between the two cell types is higher than the two homotypic ones (i.e., F LD a = F DL a > F LL a = F DD a ), cells heterogeneously mix to form a checkerboard pattern. Conversely, if the homotypic adhesions are stronger than the heterotypic ones (i.e., F LL a = F DD a > F LD a = F DL a , see the central column of Fig. 9 ), we find a spontaneous cell sorting, with the formation of small clusters of cells of the same type within the domain.
If, further, the adhesion between the light cells is larger than than the heterotypic contact interactions, which is in turn larger than the adhesion between the dark cells (i.e., F LL a > F LD a = F DL a > F DD a , see the right column of Fig. 9 ), we observe the autonomous emergence of little island of light individuals surrounded by a crew of dark cells: this phenomenon is called engulfment .
From these numerical results, it further emerges that variations in the explicit form of the repulsive part of the interaction kernels, i.e., variations in the value of the parameter s , do not significantly affect the final configuration of the system, as it is possible to see by comparing the different rows in Fig. 9 .
For the sake of completeness, we finally test the same hierarchies of intercellular adhesive/repulsive parameters, in the case of specific values no longer satisfying the H-stability criterium. As reproduced in Fig. 10 , some features of the above-observed patterns still emerge (e.g., a clusterization in little homotypic islands for F LL a = F DD a > F LD a = F DL a ): however, they are not significant from a biological point of view due to dramatic individual overlapping and aggregate collapse.
Summing up, we can therefore claim that, for any given interaction kernel (i.e., for any given value of s ), the stable minimal intercellular distance is determined by the H-stability of the system, whereas the exact cell distribution by the specific parameter setting.
These phenomena were investigated also with other types of models, see Graner and Glazier (1992) ; Mackey et al. (2014) ; Murakawa and Togashi (2015) , where the dynamics of the system are mainly driven by the hierarchy of intercellular adhesive forces, rather by their exact values.
The cell configurations obtained by our simulations are consistent with the experimental literature as well. For instance, Katsamba and coworkers in Katsamba et al. (2009) demonstrated that two identical populations of Chinese hamster ovary (CHO) cells exhibit homotypic cell sorting and form separate aggregates when express different cadherins (have different adhesive affinity). This corresponds to the numerical outcomes shown in the middle column of Fig. 9 . In other cases, i.e., when all CHO cells express the same cadherin (i.e., have the same adhesive affinity), intermixed/checkerboard aggregates are formed, in agreement with our simulations proposed in Fig. 9 , left column. Steinberg ( Steinberg and Takeichi, 1994 ) instead combined two populations of L cells transfected with P-cadherin cDNA and expressing this homophilic adhesion molecule in substantially different amounts. He then found that, when allowed to fuse, the cell population expressing more P-cadherin was enveloped by its partner, which formed an external "cortex": it is the experimental counterpart of the engulfment phenomenon observed in our results with comparable parameter setting (cf. Fig. 9 , right column). Consistently, Curtis and colleagues proposed in Curtis (1961) that the formation of internal and external cell layers during sorting results from the timing of postulated changes in cellular adhesive and motile properties.

Early migration of the zebrafish posterior lateral line (pLL) primordium
The lateral line is a sense organ, present in fish and amphibians, that is formed by a set of mechanosensory hair cells distributed in a species-specific pattern over the surface of the animal's body. These multicellular structures, named neuromasts, have the function to detect displacements and vibrations of the surrounding water: the lateral line is indeed involved in several aspects of the individual life as, for instance, prey detection, predator avoidance, and sexual courtship ( Ghysen and Dambly-Chaudière, 2004;Gompel et al., 2001 ). In zebrafish, the biological system of our interest, the lateral line extends from the head to the caudal fin and divides in two major components: the so-called anterior lateral line (aLL), which comprises the neuromasts present on the head, and the posterior lateral line (pLL), which conversely includes the neuromasts on the trunk and the tail along each side of the animal ( Ghysen and Dambly-Chaudière, 2004;Gompel et al., 2001 ). During the embryonic stage of life of the zebrafish, the pLL consists of a primordium , i.e., a proto-organ which is first recognized around 18 hpf (hours-post-fertilization), located just posterior of the otic vesicle and formed by nearly 100 epithelial cells ( Gompel et al., 2001;Haas and Gilmour, 2006 ). In particular, these component cells present several mesenchymal determinants, i.e., loss in apicobasal polarity, reduced expression of adhesive proteins and increased numbers of activated dynamic filopodia ( Lecaudey et al., 2008 ). In this respect, they are hereafter termed pseudomesenchymal individuals.
The development of the zebrafish primordium consists of two processes: (i) few hours after the formation, the proto-organ begins to migrate from the head to the tail of the embryo through the horizontal myoseptum, a connective tissue separating the dorsal and ventral body muscles masses. In particular, biological ev- Fig. 9. Final distribution of the two-population aggregate (i.e., composed of light "L" and dark "D" cells and identified by circles of radius d r ) for distinct choices both of the behavior of the interaction kernel at the origin, i.e., the value of s in Eq. (5) and of the values of adhesion strengths. All the tested parameter settings satisfy the H-stability condition. The stable pattern is mainly dependent by the hierarchy of homotypic and heterotypic cell-cell adhesiveness, whereas the value of s does not have a substantial impact.
idences show that the primordium moves along a stripe of the chemokine stromal-derived factor 1 (SDF-1). This chemical is detected only by the pseudo-mesenchymal cells located at the leading region the primordium though the expression of receptors CXCR4, while the rest of the aggregate moves upon dragging ( Donà et al., 2013;Ghysen and Dambly-Chaudière, 2004;Haas and Gilmour, 2006 ); (ii) at the same time, the activation of fibroblast growth factor (FGF) signaling pathways within the rear part of the migrating placode regulates a complex cycling process of receptor activation/differentiation that induces the organization of small groups of pseudo-mesenchymal cells into compact rosettelike structures, by increasing their epithelial character also by a loss of the migratory determinants. These clusters are then de-posited at regular intervals along the horizontal myoseptum and will mature into full active neuromasts. The rest of the primordium (still constituted by pseudo-mesenchymal cells) eventually continues its migration ( Dambly-Chaudière et al., 2007;Durdu et al., 2014;Nechiporuk and Raible, 2008 ).
Our aim is to demonstrate that the study of the characterization of H-stable intercellular interaction kernels, and of the estimate of the relative parameters, is mandatory also in the description of selected features of first of the two above-described processes, i.e., of the head-to-tail crawl of the primordium. In order to do this, let us take into account in our theoretical approach two subgroups of pseudo-mesenchymal cells: Fig. 10. Final distribution of the two-population aggregate (i.e., composed of light "L" and dark "D" cells and identified by circles of radius d r ) for distinct choices both of the behavior of the interaction kernel at the origin, i.e., the value of s in Eq. (5) and of the values of adhesion strengths. All the tested parameter settings no longer satisfy the H-stability condition. The stable pattern is mainly dependent by the hierarchy of homotypic and heterotypic cell-cell adhesiveness, whereas the value of s does not have a substantial impact.
• those that express the chemokine receptors CXCR4b, i.e., sensitive to Sdf-1. They are hereafter called leader cells (and la- sents a small portion of the zebrafish horizontal myoseptum, which is usually 30 0 0 μm long. According to the above-described experimental considerations, the general particle model described in (2) can be in this case specified as follows Fig. 11. Migration of the zebrafish pLL wild type primordium, in the case of a parameter setting that satisfies the H-stable condition for all interaction kernels K pq , i.e., for any p, q ∈ {L,T}. Cell positions are identified by colored circles with radius d r . Fig. 12. Evolution of the system with the same interaction parameters as in Fig. 11 except the ones in red, which result in disruptions of the H-stability of one or more interaction kernels. In all cases, the normal migration of the primordium is inhibited. Cell positions are identified by colored circles with radius d r .
The directional velocity characterizing the dynamics of leading cells is given by where v L is the mean migration velocity of the zebrafish protoorgan, which is approximatively 69 μm/h ( Haas and Gilmour, 2006 ), while n x is the horizontal unit vector. In this respect, we are assuming, for the sake of simplicity, that the right border of the domain identifies the tail of the animal (i.e., the target destination of the primordium) and that the myoseptum is characterized by a uniform and fixed distribution of Sdf-1, which results in a constant and persistent velocity contribution. For any pair ( pq ) ∈ {LL, LT, TL, TT}, the interaction kernels K pq introduced in (14) belong to the family of functions defined in Eq. (13) . Again, we can do some simplifications according to biological arguments: all cells forming the primordium have the same dimensions and repulsive characteristics, since they differ only for the expression of membrane receptors which have an effect only on their chemical sensitivity and adhesive properties. In this respect, we can assume F pq r = F r , s pq = s, d pq r = d r = 7 μm and d pq r = d a = 20 μm for any p, q ∈ {L, T}. In particular, cell dimensions are taken from the biological literature ( Gompel et al., 2001 ). Assuming the symmetry of F LT a = F T L a and fixed the representative value s = 1 (which implies hyperbolic repulsion at small enough inter-particle distances), we have the following set of free parameters: F r , F LL a , F LT a , F T T a ∈ R . According to the theoretical and numerical studies presented in Section 3 , realistic cell configurations can be obtained only if all the intercellu-lar interaction kernels are H-stable, i.e., if the relative parameters satisfy the criterium of Corollary 2 . In this case, given the abovedefined coefficients and hypothesis, it is indeed necessary that which results in a substantial restriction of the free parameter space. A further fitting with pLL experimental velocity and behavior allows, after a series of preliminary simulations (not shown), to opt for the following set of values F r = 1 . 12 μm / s , F LL Fig. 11 , the proposed particle model is able to reproduce a representative part of head-to-tail migration of the wild-type primordium (i.e., from the left to the right region of the domain). In particular, the leader cells guide the movement of the rest of the proto-organ only through the heterotypic adhesive interactions. In this respect, we remark that our minimal model shows good qualitative agreement with experimental outcomes without including other velocity components, successfully employed in other more specific works ( Costanzo et al., 2015 ).
To highlight the importance of the H-stability properties of the system, we then analyze the development of the computational primordium in the case of disruptions of relation (16) for one or more kernels K pq , with ( pq ) ∈ {LL,TT, LT}. This can be done either by decreasing the common particle compressibility ( F r ) or by increasing selected homotypic and/or heterotypic intercellular cell adhesiveness. The results, summarized in Fig. 12 , show that a normal collective motion of the proto-orgain is no longer obtained in the case of not completely H-stable systems. Entering in more details, we can observe that if the heterotypic interaction kernels result not H-stable (i.e., if F LT a = F T L a < F r / 27 . 21 or F LT r = F T L r < 27 . 21 F LT a ) then leader and follower cells overlap in central region of the aggregate. Conversely, if the H-stability is not satisfied from the homotypic interaction kernels, then the embryonic zebrafish pLL divides in several clusters, characterized by the collapse of particles of the same type (i.e., see the cases F LL a < F r / 27 . 21 and F LL r < 27 . 21 F LL a for the shrink of the leading region of the aggregate and the cases F T T a < F r / 27 . 21 and F T T r < 27 . 21 F T T a for the shrink of the trailing area).

Conclusions
Mathematical approaches applied to biological problems embody a wide range of techniques, which depend on the particular spatio-temporal scale of interest. However, most of them fall in two broad categories: continuous and discrete models. Continuous models approach biological phenomena in terms of variation of fields. Characteristic of a macroscopic point of view, these methods represent populations of biological individuals as densities, which evolve satisfying sets of balance laws or diffusion equations. This type of approach includes the multiphase models, developed under the simple observation that biological systems are made of several constituents and are then treated with classical concepts of continuum mechanics. On the other hand, discrete models, widely known as Individual Cell-Based Models (IBMs) or Cellular Automata (CA), approach the biological problem with a phenomenological point of view, focusing on the cell-level of abstraction and preserving the identity and the behavior of individual elements (for comprehensive reviews the reader is referred to Alber et al. (2003) ; Anderson and Rejniak (2007) ; Deutsch et al. (2007) ; Drasdo (2003b) )). Indeed, these techniques represent biological individuals, with the typical length scale of a cell, as one or a set of discrete units, with rules that describe their movements and interactions. Within this family of approaches, we can further distinguish particle-based methods, where the biological individual is represented by a material point with concentrated mass and identified by its position in space, and methods that allow for extended description of cell membrane and morphology. These latter subgroup typically relies on specific discretizations of the simulation domain with meshes can be either regular (such as square or cubic grids, as in Cellular Potts Models) or irregular (Voronoi tassellations) since each biological element is typically modelled by an ensemble of elements of the mesh.
All the above-introduced families of mathematical models have their own merits and limits. For instance, continuous techniques overlook the behavior of single individuals and also fail to describe their mutual interactions. They may therefore be unsatisfactory since, in order to deal with biological processes, it is fundamental what occurs at the scale of the single element. On the other hand, particle models are computationally efficient and easily to analyze experimentally: however, they do not allow the inclusion of intercellular chemical pathways as well as the description of morphological evolutions. Finally, methods reproducing extended and detailed cell shapes are expensive from a numerical point of view and they can only simulate several individuals at once, being inefficient in providing a general outlook of the system as a whole. In this respect, the trend of the last two decades is to create computational frameworks able to span a wide range of spatio-temporal scales with a sufficient level of accuracy, offering the advantages brought by the different methods by interface and hybridization of the different models, see Anderson and Chaplain (1998) ; Drasdo (2003a) ;Lachowicz (2002) for few examples.
In the perspective of the above-considerations, in this work we have dealt with a particle model set to reproduce the dynamics of selected cell systems. In particular, each individual has been as-sumed to move according to a first-order ODE, in the overdamped force-velocity response limit. Of the possible cell velocity components we have then focused on the term relative to direct pairwise intercellular interactions, which account for nonlocal adhesive and repulsive contributions, the former including long-range cadherin-mediated mechanisms, the latter modeling cell and nucleus resistance to compression. According to us, this is a crucial velocity contribution since the analysis of cell patterning, as well as the description of the characteristic large-time configurations of cell aggregates, is a relevant issue in developmental biology. The spatial distribution of cells mediates in fact a wide range of physiopathological phenomena, i.e., from morphogenesis to cancer invasion Ilina and Friedl, 2009 ). In this respect, the main aim of the work has been to give a consistent method that allows to facilitate the estimate of the characteristic parameters of the intercellular interaction velocity.
In more details, the intercellular interactions have been described by a proper family of pairwise interaction kernels, and relative potentials. These kernels are characterized by a repulsive part, which have different slopes according to the differentiated stiffness of cell nucleus and cytsol Ilina and Friedl, 2009 ), and by a positive parabolic trend in the attractive part. The proposed kernels are intrinsically multiparametric, being determined by a set of free coefficients relative to the extension of the interaction regions, to the intensity of the interaction forces and to the specific repulsive behavior of the particles when close enough. After some simplifications that can be made according to the problem of interest (e.g., cell dimensions), a crucial help in the analysis of the effect on the system behavior of the choice of the specific kernel (and relative coefficients) has been demonstrated to be given by the concept of H-stability. It derives from statistical mechanics and allows to predict the interparticle spacing characterizing the asymptotic configuration of a cell population. In particular, our analytical results have shown that a biologically realistic crystalline pattern of a cell aggregate can be obtained if the underlying intercellular interaction kernel, and relative potential, satisfies a proper H-stability condition. On the opposite, not H-stable kernels result in unrealistic particle collapse in our setting of cells with fixed mass.
It is however fundamental to underly that the H-stability theory is only able to reduce, for any given adhesive/repulsive kernel, the range of possible simultaneous variations of the characteristic parameters. To have a specific parameter setting, it is still necessary to perform data fitting (when experimental data are available) or proper sensitivity analysis. In this respect, we recall that the definition of H-stable kernels has implications only in the minimal distance at the equilibrium of the system particles, whose dynamics and exact distribution is entirely determined by the specific parameter values. Such considerations have been confirmed by the series of simulations proposed in this work. For instance, having the confirmation that biologically reasonable patterns of a round cell colony derive only from H-stable kernels, we have then used, in Section 3.3 , a representative experimental OVCAR spheroid to have a better coefficient setting (via data fitting).
The H-stability theory is able to assure consistent cell spacing also when dealing with more complex biological processes, involving multispecies aggregates and/or different cell velocity contributions (e.g., directional locomotion), as we have shown in Section 4 . Also in this cases, the exact cell dynamics are then determined by the specific tuning of the values of different model coefficients.
It is useful to notice that the minimal system of first-order ODEs with a given interaction potential that we have used in (4) has a well-defined limit as N → ∞ when the potential and the total mass of the system are suitably rescaled. This is usually called the mean-field limit in statistical mechanics. In fact, the mass of each point in the mean-field limit becomes negligible as N → ∞ . However, in our particle-based present approach, each agent represents a cell with a fixed mass and we are only interested in the equilibrium configurations for a fixed number of cells N . Therefore, the mean-field limit equation in our present setting of particles with fixed mass is not biologically meaningful in contrast to the topic of other works, see for instance ( Alber et al., 20 06, 20 07;Dyson et al., 2012;Lushnikov et al., 2008;Markham et al., 2013;Murray et al., 2009 ).
Our present work is based on the assumption that cell adhesive/repulsive behavior is described by the family of kernels defined in Eq. (5) . However, the analytical results hold for any other pairwise interaction functions satisfying the hypothesis of Theorem 1 . For instance, Morse-like potentials or kernels characterizing by a Gaussian profile in the repulsive part and/or by a Hooke law in the attractive one can be used. Further, except from a constant directional velocity in the last application, we have substantially focused on intercellular adhesive/repulsive interactions. This is of course an oversimplification of the biological picture. Cell migration is in fact a quite complex process involving several other mechanisms and stimuli, such as chemotaxis (i.e., cell locomotion up to gradients of a diffusible chemical field) or durotaxis (i.e., cell locomotion towards stiffer regions of the matrix environment). However, as shown in Section 4 , the inclusion of other velocity components does not change the qualitative asymptotic behavior resulting from particle interaction forces.
As main conclusion, we have demonstrated that the H-stability concept is key to choose reasonable intercellular pairwise interaction kernel and relative parameters for agent-based models with real biological implications.
In this Appendix, we provide the algebraic calculations that allow to obtain the criterium given in Corollary 2 to establish the stability properties of the family of intercellular interaction kernels introduced in Eq. (5) . In particular, we distinguish two cases according to the value of s (i.e., to the slope of K in the repulsive part nearer to the origin).
Interaction kernel with s = 1. From Eq. (7) , by setting s = 1, the interaction potential reads as: where the constants of integration C 1 , C 2 , C 3 , C 4 ∈ R are estimated in order to guarantee the continuity of the potential u ( r ) and such that lim r→ 0 u (r) = 0 . In this respect, C 4 can be taken equal to 0 and the other constants result by setting the continuity at d r /2, d r and d a . The interaction potential therefore writes as Due to Theorem 1 , the above interaction potential is H-stable when i.e., the thesis of Corollary 2 . In particular, as already explained, the values of the interaction radii d r and d a can defined according to cell phenotype: indeed, the H-stability of the system translates into a constraint on the ratio between the repulsive and adhesive interactions strengths (i.e., F r and F a ), for each value of s (for each behavior of the interaction kernel near the origin).
Interaction kernel with s = 1 (hyperbolic case) When the interaction kernel at the origin is characterized by s = 1 , the formulas can be obtained from the previous case with s = 1 by taking the limit as s → 1. This gives the following potential In this case, following the same calculation as in the previous case, we have that the H-stability Theorem 1 results in the the constraint: which is a particular case of the thesis of Corollary 2 .