A least microenvironmental uncertainty principle (LEUP) as a generative model of collective cell migration mechanisms

Collective migration is commonly observed in groups of migrating cells, in the form of swarms or aggregates. Mechanistic models have proven very useful in understanding collective cell migration. Such models, either explicitly consider the forces involved in the interaction and movement of individuals or phenomenologically define rules which mimic the observed behavior of cells. However, mechanisms leading to collective migration are varied and specific to the type of cells involved. Additionally, the precise and complete dynamics of many important chemomechanical factors influencing cell movement, from signalling pathways to substrate sensing, are typically either too complex or largely unknown. The question is how to make quantitative/qualitative predictions of collective behavior without exact mechanistic knowledge. Here we propose the least microenvironmental uncertainty principle (LEUP) that may serve as a generative model of collective migration without precise incorporation of full mechanistic details. Using statistical physics tools, we show that the famous Vicsek model is a special case of LEUP. Finally, to test the biological applicability of our theory, we apply LEUP to construct a model of the collective behavior of spherical Serratia marcescens bacteria, where the underlying migration mechanisms remain elusive.

. Applicability of different modelling/computational approaches according to their biological interpretability and the corresponding knowledge of biophysical migration mechanisms. The latter is typically encoded as the degree of interaction potential knowledge in Langevin equations (see text). In the case of extended biophysical mechanism knowledge, mechanistic models are the natural choice. When the effects of cell-cell interaction on cell migration are only partially understood, phenomenological models can be typically used. Finally, when data do not suffice to formulate an interaction potential, machine learning allows for the quantitative reproduction of experimental data. However, this has a toll in the interpretability of the resulting model, since machine learning methods are typically "black boxes". LEUP models offer a compromise that allows for quantitative predictions under lack of mechanism knowledge and satisfactory biological interpretability.
Scientific Reports | (2020) 10:22371 | https://doi.org/10.1038/s41598-020-79119-y www.nature.com/scientificreports/ In this work, we present the simplest LEUP-driven Langevin model of swarming where individuals can sense the velocity orientations of other individuals in their surroundings. Individuals act as Bayesian inferrers and change their own orientation to optimize their prior, according to environmental orientation information. Under these assumptions, individuals reorient according to the entropy gradient of the environmental information. A parameter, named the sensitivity, controls the strength and directionality of the reorientation in relation to the local gradient. We find that the system adopts a steady, polar-ordered state for negative values of the sensitivity. Conversely, the system remains out of equilibrium, but partially nematic-ordered when the sensitivity is positive. Furthermore, we find that the qualitative behavior of the model depends on the values of the particle density, noise strength, sensitivity, and size of the interaction neighborhood. Finally, we showcase the LEUP principle by showing that our model replicates the collective behavior of spherical S. marcescens bacteria.

Materials and methods
The self-propelled particle framework. Moving and interacting cells are modeled by a two-dimensional self-propelled particle model (SPP). In this model, N ∈ N cells move on a two-dimensional area. The n-th cell is characterized by its position, � r n ∈ R 2 , speed, v n ∈ [0, ∞) ⊂ R , and an orientation θ n ∈ [0, 2π) ⊂ R . Due to the small size of cells, it is assumed that viscous forces dominate. Changes in speed and orientation result from local potentials U θ (� r m , θ m , v m ), U v (� r m , θ m , v m ) : R 2 × [0, 2π) × [0, ∞) � → R which depend on the positions and polar velocity components of cells within a radius R ∈ R + . The bias of the cell to follow the potential gradients are regulated by the parameters β θ , β v ∈ R , called angular and radial sensitivities, respectively. Additionally, velocity fluctuations occur due stochastic noise terms ξ α n (t) ∈ [0, 2π) , α ∈ {θ , v} where t ∈ R + denotes time. The noise will be assumed to be a zero-mean, white noise term, which has the statistical properties �ξ α n (t)� = 0 and �ξ α n (t 1 )ξ α m (t 2 )� = 2D α δ(t 1 − t 2 )δ nm , where t 1 and t 2 are two time points, D α ∈ R + is either the angular ( α = θ ) or radial ( α = v ) diffusion coefficient, δ(t) is the Dirac delta, and δ nm is the Kronecker delta. Finally, the radial acceleration will be assumed to be damped by a density dependent friction, ψ(ρ n ) . In the following, it will be assumed that the density-dependent friction is given by ψ(ρ n ) = ρ n −ρ , where ρ n is the local cell density within the n-th cell's interaction radius, and ρ is the global average cell density. Taking everything into account, the stochastic equations of motion of the n-th cell read 16 where v(θ n ) is the normalized velocity of the cell and ε is a parameter. A representation of the SPP model is shown in Fig. 2. The function g( v n ) modulates the noise variance and allows us to model certain distributions (such as the Rayleigh distribution in Section 5).

Figure 2.
Graphic representation of the dynamics of the SPP model. The n-th cell is represented by a point particle with speed v n and orientation θ n . Depending on the form of the interaction potential, the cell may feel a reorientation force −∂ θ U θ and a radial force −∂ v U v due to interaction with other cells inside the interaction neighborhood defined by the radius R.
Scientific Reports | (2020) 10:22371 | https://doi.org/10.1038/s41598-020-79119-y www.nature.com/scientificreports/ The interaction potentials U α ( r m , θ m , v m ) , which dictate the velocity dynamics of cells, need to be specified. Biophysically, the potentials should encompass steric effects, hydrodynamic interactions, chemotactic effects, and terms arising from internal cellular processes, for example, flagellar motor dynamics, actin polymerization, receptor dynamics, etc. Finding such potentials is a formidable task since not all of the mechanisms and interactions involved are known. To circumvent this problem, a variational principle of cell decision-making related to entropy maximization 17 , known as the least microenvironmental uncertainty principle (LEUP), will be used 18 . In the next section, we will discuss such a case.
Least microenvironmental uncertainty principle (LEUP). The main premise of LEUP is that cells equip Bayesian inference to decide their internal states, expressed as a combination of a sensed microenvironmental distribution (empirical/measured likelihood distribution) and a local entropy-dependent prior. From the cellular point of view, constructing an accurate microenvironmental sensing distribution is expensive, since cell sensing is an energetically costly process. Cells sense (collect information) of their surroundings by employing different processes, such as polymerizing pseudopodia, translocating receptor molecules or modifying its cytoskeleton according to mechanical signals 19,20 . However, the cost can be minimized when cells build informative priors about their microenvironment. In the case of moving cells this could be achieved by, for example, promoting cell polarization through the recruitment of actin related proteins 21 . In turn, this allows cells to spare the energy from building extra sensory processes. At same time, cell polarization effectively promotes the coevolution of the cellular state and its microenvironment, where the latter becomes more predictable. The simultaneous determination of cell state and its microenvironment results in the minimization of the microenvironmental uncertainty towards a target value related to a given tissue or bacteria (e.g. biofilm). This process also is prominent during differentiation, where pluripotent progenitors generically sense, with the goal to find an appropriate differentiation niche, but differentiated cells have very targeted sensors, e.g. precise receptors, that detect a narrow range of tissue-relevant cues. This gives them an advantage since they can invest their energy in optimizing their actual function rather than environmental sensing. The state of the n-th cell in this case is defined by its orientation θ n and velocity v n . We assume that the orientation and velocity of cells are decoupled, i.e. one can consider orientations and velocities independent from one another. The set of intrinsic angular states of other cells within its radius of interaction is given by � n = {θ m : �� r n − � r m � ≤ R} , while the set of intrinsic velocity states is V n = {v m : �� r n − � r m � ≤ R} . The cell reacts to the environmental information, n and V n , by changing its own states, θ n and v n . The cell then acts as a Bayesian decision-maker, such that α n ∈ {θ n , v n } , A n ∈ { n , V n } where P(A n | α n ) can be interpreted as the accuracy with which a cell can sense other cells in their surroundings and react accordingly, and P(α n ) is the probability distribution of the cell's intrinsic states (or prior). However, sensing other cells and evaluating P(A n | α n ) entails an energy cost. It is reasonable to assume that the cell will try to optimize its prior P(α n ) for the sake of energetic frugality.
The prior probabilities P(α n ) should, consequently, fulfill certain the premises of LEUP. First and foremost, they should be normalized, i.e. P(α n = a)da = 1 , integrating over all possible values of internal cellular states. Second, a biological cell is an imperfect sensor. Therefore, the uncertainty in sensing accuracy S(A n | α n ) , should reach a certain level in average, which is species-dependent. This assumption is a cornerstone of the LEUP formalism. An important feature of the LEUP formalism, and of the entropy maximization principle in general, is its applicability when every other mechanistic details of the system are unknown. Entropy is a measure of uncertainty according to information theory, therefore entropy should be maximal, in order to reflect our lack of mechanistic knowledge of the phenomenon, and to avoid introducing any artificial bias in the model arising from the specific choice of P(α n ) 22 . Given that the internal entropy is given by S(α n ) = − P(α n = a) ln P(α n = a)da , entropy maximization subjected to probability normalization and the target mean sensing accuracy translates into the optimization problem where δ δP(α n ) is the functional derivative, S (A n | α n ) is the target sensing accuracy, and and β α are Lagrange multipliers. Taking into account the relations among entropy and probability, Eq. (2) yields where Z = e −β α S(A n |α n =a) da is a normalization constant and β α is the responsiveness of the cell. Using Eq. (3), the internal entropy of the cell, defined as S(α n ) = − P(α n = a) ln P(α n = a)da , is given by Using the relation between thermodynamic-like potentials, it is evident that the average internal energy is given by Scientific Reports | (2020) 10:22371 | https://doi.org/10.1038/s41598-020-79119-y www.nature.com/scientificreports/ and for a particular realisation of α n we can write the internal energy as Now, Helmholtz-like free energy can be written as Finally, we return back to the Bayesian formalism and to the main LEUP premise that cells tend to build informative microenvironmental priors P(α n ) . In the steady state, the latter implies that i.e. the posterior distribution of P(α n |A n ) becomes independent of measuring/sensing the microenvironmental information A n , since the prior P(α n ) includes all relevant information. Observing the Bayesian probability, this only happens when the measurement/sensing likelihood P(A n |α n ) approaches the probability distribution of the microenvironment P(A n ) , i.e. the cell perfectly senses its microenvironment. Using again the Bayesian probability, we can prove the following information-theoretic equation Substituting the above in Eq. (4), the mutual information reads Two cases have to be considered, depending on the sign of β α : • When β α < 0 the only positive term is ln Z in Eq. (10). Thus, it is reasonable to expect the mutual information to decrease to a minimum value. • When β α ≥ 0 one can show that the system, for a large enough sensing radius, goes to an ordered state at the steady state (see "Collective cell migration patterns for different parameter regimes"). This implies that the microenvironmental entropy S(A n |α n ) tends to zero. In this case, the Eq. (10) becomes This implies that the mutual information has the potential to reach zero. Interestingly, if this happens using Eq. (9), we can recover a similar definition to Boltzmann entropy for the internal cell state entropy, i.e.
LEUP-based dynamics. The internal energy depends on the internal states of the cell, as well as the internal states of other cells in the surroundings. Such internal states can be a vector of physical quantities (e.g. velocity, acceleration) and/or chemical variables such as intracellular proteins, genes and so on. Here, we focus on the former to define an interaction potential that models the equations of motion. By doing so, it is evident that the responsiveness of the cells to LEUP can be quantified by the sensitivity β α = −β α . Analogously, we can write the equations of motion of the model similar to Eq. (1c) as To illustrate entropy calculation, it will be assumed that the orientations of cells within the interaction neighborhood are distributed according to where µ is the mean of the distribution and γ is a parameter related to the variance. This is a wrapped Cauchy distribution, periodic over the interval [0, 2π] . Similarly, cell speeds will be assumed to be distributed half-normally www.nature.com/scientificreports/ where σ 2 is proportional to the variance of the distribution. Accordingly, the angular entropy is while the speed entropy is The parameter σ can be determined from the local speed variance, while the parameter γ depends on the local polar order (i.e. the degree of parallel alignment) of cell velocities in the neighborhood. It should be noted that the qualitative behavior of the model is independent of the particular choice of distributions, and the distributions considered here are suggested only for ease of calculation. Before defining γ , we will first define the observables characterizing the order of the velocity field.
Collective migration observables. Let us define the normalized complex velocity of the n-th cell, z n ∈ C as z n = e iθ n , where i is the imaginary unit. The k-th moment of the velocity over an area A is given by where the sum is over all cells in area A, and N A is the total number of cells in A. The polar order parameter in the area A is given by which is the modulus of the first moment of the complex velocity in A, while the nematic order parameter in the area A is given by which is the modulus of the second moment of the complex velocity in A. The order parameters are bounded, i.e.
due to the complex velocities z n being normalized. The parameter γ for the distribution of orientations in the neighborhood of the n-th cell is given by where the subindices C R,n indicate a circular area of radius R centered at r n . The latter directly stems from the properties of the wrapped Cauchy distribution.
While global polar and/or nematic order are characteristic of steady flows, rotating flow fields are commonly observed in out-of-equilibrium systems. The vorticity is an observable which is equal to twice the local angular velocity, and is thus a measure of the local strength and direction of rotation of the field. The vorticity ω is defined as where v mean ( r) is the mean velocity field at point r , and k is vector normal to the plane where cells move.

Statistical evaluation of experiments and model predictions.
can be viewed as a z-score for each Ô (j) i and for large enough simulation ensemble should converge to a normal distribution. The total degrees of freedom for both observables is 2N = 22 . The we calculate the reduced χ 2 −statistic, or χ 2 per degree of freedom, which is defined as χ 2 2N = (2N) −1 j χ 2 j = 1.97 being close enough to 1. This suggests that our fitting is satisfactory, since values χ 2 2N ≫ 1 indicate a bad fit to the experimental data.

Results
From LEUP to phenomenological models of collective migration: the relationship to the Vicsek model. By using the LEUP, we have modeled interaction as a change in velocity dictated by the local entropy gradient. The modulation of β α parameters modulates the response of cells to the local entropy gradient and gives rise to relationships with known phenomenological models, such as the Vicsek model. The absolute value |β α | is proportional to the likelihood of the cell to change its velocity according to a given entropy gradient. If β α < 0 , cells tend to go against the local entropy gradient towards the entropy minimum. In the specific (14) S(� n | θ n ) = ln(2π) + ln(1 − e −2γ ), Scientific Reports | (2020) 10:22371 | https://doi.org/10.1038/s41598-020-79119-y www.nature.com/scientificreports/ case of α = θ , a negative sensitivity would restrict the distribution of angles to a narrow selection. Conversely, β α > 0 forces cells to follow the entropy gradient towards the entropy maximum, broadening the distribution. From here on, we will assume that the effect of cell interactions will be averaging the radial component, therefore To evaluate the effect of these two opposite migration strategies, we analyze the angular steady states in the two parameter regimes. Without loss of generality, we assume that, in the steady state, the mean velocity is v = 1 . By expanding S 1 C R,n using Eq. (16), defining the components of the mean neighborhood velocity as v y,n = C R,n ∋m� =n sin θ m and v x,n = C R,n ∋m� =n cos θ m , and differentiating Eq. (14), we find that the orientation of θ n at the entropy extrema must be such that (see Supporting Information) but v y,n v x,n = tanθ , the tangent of the mean orientation of the neighbors, excluding the n-th cell. This results in two extremum points θ n =θ and θ n =θ + π , one where the velocity of the n-th cell is parallel to the average velocity of its neighbors, and one when it is antiparallel. In the first case while in the second case It can be shown (see Supporting information) that θ n =θ corresponds to an entropy minimum, while θ n =θ + π corresponds to an entropy maximum. Consequently, the behavior of the regime β θ < 0 is analogous to that of the Vicsek model 6 . Conversely, the regime β θ > 0 corresponds to an anti-ferromagnetic analog of the Vicsek model.
Next, let us assume that the model has a steady state, where the Helmholtz free energy per cell is given by Eq. (7). Due to its extensivity, the Helmholtz free energy of complete, non-interacting, steady state system is where Z n is the normalization constant of Eq. (3) for the n-th cell. For a weakly interacting system, the mean-field effective normalization constant Z T := N n=1 Z n is given by Note that this is only valid in the limit β θ → 0 . Integrating and substituting the resulting Z T into Eq. (7) (see Supporting information), yields the Helmoltz-like free energy Eq. (25) is well-defined only for β θ < 1 . This indicates that no steady state exists for β θ ≥ 1 , hinting at an out-ofequilibrium regime 23 . The present model belongs to the class of models with logarithmic potentials (see Eqs. (6) and (14)). The existence of a non-normalizable state in certain parameter regimes is a staple of systems with logarithmic potentials 24 .

Collective cell migration patterns for different parameter regimes. The model was implemented
computationally to characterize the model and the effects of the different parameters on the resulting macroscopic behavior. The general qualitative behavior of the model can be observed in Fig. 3. In the regime β θ < 0 , cells tend to travel in a single direction after some time has elapsed, similar to the Vicsek model. Conversely, in the β θ > 0 regime, cells are seen to move collectively in transient vortex-like structures, even after long times have elapsed. Qualitatively, the patterns resulting from different parameter combinations are summarized in Table 1. Analyzing simulations, two important phenomena are observed. First, there is a critical parametric regime � C := {(β θ , R) : S 1 A , S 2 A > 0} where patterns emerge. Specifically, for low values of interaction radius R no structures can be formed. This indicates that medium-to-long range spread of information is necessary for ordering. On the other hand, for β θ > 0 for large values of R, outside of C , again no patterns occur. This implies that for large interaction radii there is a destructive interference of the travelling information. A second important observation is that patterns do not depend on the choice of β v when this is different than zero. For β v = 0 , LEUP dynamics divide the population into fast and slow cells. While fast cells are useful for spreading information (and therefore, increasing the effective interaction range), slow cells are necessary for maintaining local ordering. On the other hand, if we fix the initial speed distribution and assume β v = 0 , then we find different patterns emerging as shown in Table 1 (also see in SI). Furthermore, we quantitatively characterized global ordering at long times. The global polar order parameter, given by Eq. (16), for the complete simulation domain, tan θ n =v y,n v x,n , (22a) sin θ n ∝v y,n and (22b) cos θ n ∝v x,n ,  Table 1. Qualitative description of the observed patterns for different angular sensitivity and interaction radius regimes, as well as radial sensitivity.The patterning regime C is the blue area in Fig. 4a,b.

Polar aligned streets of cells Scattered polar aligned cells
β v = 0 (for uniform distribution) Compact polar aligned cluster Compact polar aligned cluster No order or patterns Vortices   17) for the complete simulation domain, measures the tendency of all cells to align nematically, or along a single axis. These order parameters take a value of one when there is global order, while taking a value of zero when the system is completely disordered. It should be noted that polar order implies nematic order, but the reverse is not true. Similarly to other velocity alignment models 7 , the model shows an order-disorder transition with increasing noise amplitude and decreasing density (SI figure Fig. 1). More importantly, we observe that in the regime β θ < 0 the system also undergoes a transition towards polar order with decreasing β θ . After the transition, most particles have a similar orientation (Fig. 5a,c). In the regime β θ > 0 , a phase transition is also observed towards nematic order with increasing β θ . In this case, however, the nematic ordering is not perfect, as evidenced by the nematic order parameter reaching values of around 0.35 after transition (b) compared to the value of 0.9 of the polar order parameter after transition in the β θ < 0 regime. This is further evidenced by the bimodal distribution of orientations with peak separation of approximately π radians (Fig. 5d). These simulation results further corroborate our previous theoretical results.
In turn, we study the effect of speed sensitivity β v in terms of phase transitions. We fix the angular sensitivity β θ , either positive or negative, thus the speed distribution will only depend on β v values.
When β v is positive then speed distribution become bimodal and for β v < 0 the speed distribution becomes unimodal (see Fig. 10 in supplementary material). Moreover, we show that if the radial sensitivity β v < 0 decreases then the average speed increases. On the other hand, for any value of β v > 0 , the first and the second www.nature.com/scientificreports/ moments of the speed distribution cannot be defined, since this is bimodal. Finally, for increasing cell densities, the average speed increases as well (see supplementary Fig. 9).
A collective migration example of restricted mechanistic knowledge: the spherical bacteria case. Collective motion of bacteria has been extensively studied and modeled. Most studies have focused on the collective properties of S. enterica, E. coli, and M. xanthus. These species of bacteria are similar since they have a high aspect ratio. It has been shown that volume exclusion, coupled with a high aspect ratio, is sufficient to induce velocity alignment in the system 7 , and accordingly, ordered clusters of bacteria are observed at high densities. However, it has been recently shown 25 that even spherical S. marcescens bacteria do display collective migration (for experimental details please see SI section). The biophysical mechanism whereby spherical bacteria interact with one another must be different from the high body aspect ratio volume exclusion mechanism proposed for elongated bacterial species.
Recently, a combination of biophysical agent-based and hydrodynamics model has been proposed to describe these experiments. In this study the experimental observations were only partially reproduced. Therefore, the biophysical mechanisms underlying collective migration in spherical bacteria are still not well understood. An important aspect to consider is the bacterial speed v n . It was found experimentally 25 that bacterial speed followed a Rayleigh distribution, dependent on bacterial density. Collective effects on cell orientations, on the other hand, were studied by observing the vortical behavior of the population 25 .
To reproduce the experimentally observed Rayleigh distribution for cell speed, we chose the function g( v n ) = v −1 n as shown in 16 . It is important to note that this term is not impacting the qualitative behavior of the average bacteria speed but only its variance (see SI Fig. 11). Moreover, the interested reader could see the impact of the friction term in the average cell velocity in SI Fig. 12. As shown in Fig. 6, our model qualitatively and quantitatively reproduces both the speed distribution and vorticity behavior of the experimental system. Interestingly, the behavior of the experimental system was replicated for high values of the sensitivities β v and β θ , and large interaction radii R.
Our LEUP model not only allows for a quantitative reproduction of the experiments, but also provides insight into the potential biophysical mechanisms. Such values of the sensitivities and interaction radii indicate far-reaching, strong tendencies of bacteria to average their speeds while reorienting and traveling differently from their neighbors. Spherical, rear-propelled particles have been shown to destroy polar order as a result of hydrodynamic interactions 26 , similarly to our model. Considering that S. marcescens is an example of a spherical, rear-propelled particle 27 , our results agree with previous findings indicating that S. marcescens interacts through long-range hydrodynamics 28 . The long range interaction radius suggests the existence of hydrodynamically induced interaction (which has been suggested by Ariel et. al as well as by other studies 27,28 ) or self avoiding interaction 29 .

Discussion
In this work, we have introduced an off-lattice model of LEUP-induced collective migration, based on the self-propelled particles modeling framework. It was assumed that individuals changed their radial and angular velocity components independently through LEUP. Reorientation is governed by a stochastic differential equation depending on a white noise term and a force arising from an interaction potential.
The exact form of the interaction potential can be very complex, and its specific form is dependent on particular mechanochemical details of the modeled system. While it has been shown that, in general, interactions among individuals can effectively drive the entropy of the entire system towards an extremum point 30,31 , here we do the opposite. Instead of modeling the interaction potential biophysically, it was assumed that particles followed the LEUP, which dictates that cells change their internal states in order to minimize the uncertainty of the internal states of cells in their surroundings. Although LEUP has been conceptualized to deal with highdimensional internal states involved in cell decision-making, here we restrict on physical internal states such as speed and orientation. While cell speed was assumed to always minimize uncertainty, there was no assumption made on the cell orientation. Particles are therefore free to reorient either towards or against the gradient of entropy of the orientational distribution of particles in their neighborhood, depending on the sign of the sensitivity parameter, which also dictates the strength of the interaction. The orientational distribution in the neighborhood was assumed to be wrapped Cauchy distributed. Such a distribution facilitates the mathematical analysis of the model. However, the usage of other wrapped distributions do not qualitatively change the general behavior of the model (see SI). Please note that non-parametric methods for estimating entropies without assuming any underlying parametric distributions exist. For instance, such methods employ kernel density estimation, k−nearest neighbours or regression methods 32 .
We show that, when the parameter β θ is negative, the model produces steady-state polar alignment patterns. Interestingly, we showed that the classical formulation of Vicsek model 6 is a special case of LEUP. Conversely, when the parameter β θ is positive, particles tend to reorient against the mean velocity of their neighborhood. In this regime, the free energy diverges, indicating an out-of-equilibrium parameter regime. This kind of parameter-dependent dichotomy is similarly observed in systems with logarithmic potentials 33 , involved in processes such as long range-interacting gases 34 , optical lattices 35 , and DNA denaturation 36 . The dichotomy arises from the logarithmic form of the entropy driving interaction in our model. It has been shown that, due to the nonnormalizability of the steady state solution, such systems require a time-dependent expression for their analysis 24 . Therefore, an in-depth theoretical analysis of our model would require a similar multiparticle, time-dependent expression of the angular probability densities.  37 we have developed a discrete speed version of our LEUP migration model, where cells can have only zero or a finite speed. This model exhibits Turing patterns, i.e. dynamics clusters of non-motile cells of specific characteristic wavelength, where previously published Viscek-like models cannot may produce moving clusters of swirling cells (e.g. the milling Viscek model) but never static ones.
As a proof of principle, we show that our model replicates the collective vortical behavior of spherical motile particles. Recently, the collective behavior of spherical particles have been modeled as a combination of steric repulsion and hydrodynamic interactions 38 . Our study has shown that hydrodynamics and steric interactions induce long-range microenvironmental entropy maximization, which coincides with the β θ > 0 LEUP regime. This generalizes the type of biophysical mechanisms required to produce vortical patterns.
It should be noted that, while spherical S. marcescens bacteria have been modeled biophysically, their collective behavior was partially reproduced 8 . This hints at an additional biological and/or biochemical interaction between cells. While our LEUP-based model is coarse-grained in terms of specific biophysical/biochemical interactions, it allows for a plausible reproduction of the experimentally observed collective velocity behavior by fitting a only few parameters. The application to spherical bacteria allows us to showcase the potential of the LEUP principle when the precise interaction mechanisms are not known.
As already mentioned, we have made some assumptions to simplify the model. Our model assumes a Gaussian, white noise term in the SDEs. This results in normal diffusive behavior in the absence of interactions. It has been observed experimentally, however, that in some conditions, cells perform Lévy walks resulting in  25 . Throughout all simulations, the standard deviation of the noise was set at 0.0001, interaction radius at R = 10 , proportionality constant ε = 0.008 , radial sensitivity β v = −20 , g = 1 v n and angular sensitivity at β θ = 5 . Data was obtained after 500 time steps.

Scientific Reports
| (2020) 10:22371 | https://doi.org/10.1038/s41598-020-79119-y www.nature.com/scientificreports/ superdiffusive behavior 39 . By changing the distribution or time correlations of the noise 9,40 , it would be possible to both replicate the non-Gaussian dynamics of single cells, and investigate the effect of single anomalous dynamics on collective behavior. We have also assumed that particle velocities are the only internal states relevant for reorientation, for simplicity and as a proof of concept of the LEUP principle. However, it is reasonable to think that other states, such as relative position or adhesive state, may be relevant to include when modeling specific systems. This reveals an interesting point in the application of LEUP-driven models which is the selection of the most relevant/dominant internal variables. Although experimental intuition could be the easiest approach, we are currently developing a spatial principle component analysis method that would allow to select the most relevant internal variables using spatial data such as multiplexing biopsies or spatial RNA sequencing.
As stated above, LEUP circumvents the biophysical details of cell migration. The need to model systems of interacting agents without previous knowledge of the biophysical mechanisms involved has sparked at least another agent based model 41 . In this model, similarly to ours, agents act without a mechanistic rule. Rather, they consider every possible action and penalize those which are not favorable to their internal standards. While both the aforementioned model and LEUP are defined in a similar spirit, modeling under LEUP consists in correctly identifying the relevant internal cellular states for entropy optimization, while in 41 modeling is concerned with defining suitable penalizations for each possible decision scenario.
LEUP has additional appealing features. For instance, LEUP allows for replicating a plethora of collective migration patterns. In this particular case, we have analytically derived the polar and nematic alignment Vicsek models for LEUP arguments. In this sense, LEUP acts as a generative model for collective migration mechanisms. This is particularly useful upon limited knowledge of such mechanisms, a problem called structural model uncertainty. Another advantage of LEUP is the mapping of biophysical mechanism combination to the β > 0 or β < 0 regimes. This allows for unifying the model analysis but for a better classification of migration mechanisms. Finally, known mechanisms or data could be easily integrated to our proposed framework by further constraining the LEUP dynamics.