Simona Mancini
email: simona.mancini@univ-orleans.fr

René-Marc Mège

Benoit Sarels

Pierre-Olivier Strale

R.-M Mège

P.-O Strale

A phenomenological model of cell-cell adhesion mediated by cadherins

Keywords: Cell Adhesion, Cadherins, Adherens Junctions, Protein Clustering, Degenerate Reaction-Diffusion System, Patterns Formation Mathematics Subject Classification (2000) MSC 92C15, MSC 82B21, MSC 35J70

We present a phenomenological model intended to describe at the protein population level the formation of cell-cell junctions by the local recruitment of homophilic cadherin adhesion receptors. This modeling may have a much wider implication in biological processes since many adhesion receptors, channel proteins and other membrane-born proteins associate in clusters or oligomers at the cell surface. Mathematically, it consists in a degenerate reaction-diffusion system of two partial differential equations modeling the time-space evolution of two cadherin populations over a surface: the first one represents the diffusing cadherins and the second one concerns the fixed ones. After discussing the stability of the solutions of the model, we perform numerical simulations and show relevant analogies with experimental results. In particular, we show patterns or aggregates formation for a certain set of parameters. Moreover, perturbing the stationary solution, both density populations converge in large times to some saturation level. Finally, an exponential rate of convergence is numerically obtained and is shown to be in agreement, for a suitable set of parameters, with the one obtained in some in vitro experiments.

Introduction

Intercellular junctions are macromolecular structures built at the interface between cell membranes and holding animal cells together within a tissue. Two of these types of junctions, adherens junctions and desmosomes, are formed by the local recruitment of transmembrane proteins of the cadherin family tightly associated intracellularly to the cytoskeleton which provides cells with particular viscoelastic properties and mechanical resistance. Adherens junction formation deserves much interest among experimental biologists, since it initiates the formation of all the other types of intercellular junctions, including desmosomes, and thus it is at the centre of the cohesion and mechanical resistance of biological tissues. Cell-cell adhesion associated to adherens junction formation is initiated by the homophilic trans-interaction of extracellular domain of cadherins from two adjacent cells, see [?]. The strengthening of adhesion is then favored by the clustering of cadherin molecules in the two membranes and by the anchoring of their cytoplasmic domain, via catenin adapter proteins, to the underlying actin filaments (F-actin), themselves maintained under tension by non-muscle myosin motors (myosin II). However, because of the high variability of cell-cell contact shapes and because of the multitude of partners and regulatory steps involved, the sequence of molecular and cellular events leading to the formation of mature adherens junctions from this cellcell contact initiation step is still unclear. Thus, it is important to understand how the trans-interaction of cadherins leads to their clustering and anchoring to F-actin, so to form discrete size-defined junctions, see [?,?,?]. To overcome these difficulties we described earlier an experimental approach where single cells are allowed to spread on surfaces covered with purified cadherin extracellular domains, mimicking the cell-cell contact formation, see [?]. On such a substratum, cells adopt an isotropic morphology with radial accumulations of cadherin adhesion complexes on their ventral face that colocalize with F-actin, showing the formation of adherens junctions (Fig. 1).

Interestingly, these local accumulations of cadherins also form when the F-actin network is destructured by pharmacological means, indicating that cadherin clusters formation could be dictated solely by trans-interactions of cadherin ectodomains, see [?]. However, how such interactions could generate the observed patterned distribution of cadherins is not obvious. Interestingly, cadherin dense regions correspond to areas of strong membrane apposition to the adhesive surface, see [?]. In between cadherin accumulations, the membrane is significantly remote from the substratum, likely as a result of either the repulsion exerted by larger extracellular glycoproteins or of spontaneously occurring plasma membrane fluctuation, see [?,?]. Further, it has been experimentally observed in [?] that cadherins uniformly bind over the whole domain as soon as the unbinding is weak.

In the present study we propose a mathematical model describing the contribution of cadherin trans-interactions to their local recruitment in adhesion aggregates. Our goal is to recover the experimentally observed behavior applying a rather simple mathematical description. In particular, our main aim is to show how, from an uniform distribution of diffusing cadherins and some bonds which number and density are negligible, we can (or not) obtain cadherin aggregates. Further, we want to analytically recover the exponential growth of the proportion of fixed cadherins (i.e. the Immobile Fraction) towards a saturation level, and its speed. The model takes into account the diffusive properties of cadherins and their binding and unbinding probabilities. It is a macroscopic model describing the time and space evolution of two density distribution functions, one for the freely diffusing cadherins and one for the fixed ones. Since one of the two population of cadherins does not diffuse and since the two populations interact, the system of partial differential equations we consider is a degenerate reaction-diffusion system. It will be numerically shown that the solution of the system is strongly dependent on the initial density distributions, and we shall therefore validate our model comparing, rather than the distribution of aggregates, the time evolution of the proportion of free and bound cadherins, i.e. the Mobile and Immobile Fractions, with those obtained experimentally.

The paper is organized as follows. In Sect. 2 we summarize the basic biological features and define the mathematical model. Its mathematical analysis is performed in Sect. 3: after proving that the stationary solution of the considered system are homogeneous in space, we perform a stability analysis to prove that taking the parameters in some range leads to structure formations. Sect. 4 is devoted to the numerical approach: first we show the behavior of the solution on some theoretical tests, then we consider more biological relevant tests and compare the numerical results with experimental ones. Some discussions are finally given in Sect. 5.

The mathematical model

We aim at modeling the clustering of cadherins seen at the macroscopic scale on cells spread on cadherin-coated surfaces (see Fig. 1). This problem is thus bi-dimensional and we consider a bounded domain Ω ⊂ R 2 . In a first approximation we simplify the problem making the following assumptions :

-On the substrate, there is a coating of cadherins that mimics the surface of a neighboring cell. We shall call these cadherins the targets and we assume that their number is constant in time and that they are uniformly distributed on the domain Ω. We describe them in the mathematical model by their density ρ > 0. Note that these targets are fixed, hence they do not diffuse.

-Over the substrate is the cell membrane on which cadherins diffuse and may bind to the substrate. The distance between the substrate and the membrane being negligible, we may say that the considered cadherins evolve on the substrate itself. We consider two populations of cadherins : those diffusing, whose density distribution will be denoted by u = u(x, t) and which we refer to as free or diffusing cadherins ; and those linked to the substrate, whose density distribution will be denoted by v(x, t) and which we will call fixed cadherins. We note that, from a biological point of view, it is not possible to have more fixed cadherins than targets, so that the density v(x, t) must be smaller than the density ρ at each time t and each position x. -A diffusing cadherin can bind to a target cadherin via trans-interactions.

The probability for a diffusing cadherin to bind is locally increased by the presence of other fixed cadherins, see [?,?]. This will be represented by a non-linear increasing, positive and bounded function F (v), which we will refer to as the binding function. This is not the only possible definition for the binding function. In fact, we may define a binding function considering that the cell membrane fluctuates over the surface, so that the higher the distance h of the membrane to the surface is, the smaller will be the probability of cadherins to bind. How to define such a distance h without introducing too many variables in the model ? We overcome this difficulty by noting that the larger is h, the less fixed cadherins can influence the behavior of diffusing one, and thus by defining the binding function F as a function of v and considering that cadherins diffuse on the same surface of targets. -Conversely, fixed cadherins may unbind. The probability for a fixed cadherin to unbind locally depends on the presence of other fixed cadherins:

the more bonds are locally present, the harder the cadherin unbind. In fact, local bonds stabilize the aggregates, for example by means of the cis-contacts. This probability will be represented by a non-linear decreasing, positive and bounded function G(v), denoted also as the unbinding function.

-As it is shown in [?], once cadherins are bound to targets, their diffusion coefficient loses two orders of magnitude. We then assume that free cadherins diffuse on the membrane while fixed cadherins do not.

Thus, the mathematical model describing the biological phenomena we consider can be described, at a macroscopic level, by the following degenerate reaction-diffusion system,

∂ t u = σ ∆u -εr(u, v) ∂ t v = εr(u, v) (1)
where σ is the diffusion coefficient and ε is a parameter describing the efficacy of the reaction term on the evolution in terms of a frequency. System (1) is composed of a reaction-diffusion equation and a reaction equation (a simple time evolution partial differential equation), with unknowns the densities u = u(x, t) and v = v(x, t), respectively representing the distributions of free and bound cadherins at time t ≥ 0 and in a position

x ∈ Ω. Since Ω ⊂ R 2 ,
x is a two-component vector, (x 1 , x 2 ), and ∆ is the Laplacian operator with respect to the two components of x. Note that u and v being two densities, they must be non-negative at all time, that is u(x, t) ≥ 0 and v(x, t) ≥ 0 for all t ≥ 0 and x ∈ Ω. Moreover, as previously explained, from a biological point of view, it is not possible to have more fixed cadherins than the possible available targets, thus it must be v(x, t) ≤ ρ for all t ≥ 0 and x ∈ Ω, too. Note that no upper bound is required for u(x, t).

By the first equation in (1), we have that free cadherins diffuse and may change their status, becoming fixed, following the reaction term -r(u, v). Symmetrically, the second equation in (1) says that fixed cadherins may change their status becoming free following the reaction term r(u, v). Hence, r(u, v) must describe the gain, r + (u, v), and loss, r -(u, v), rates of fixed particles.

The gain rate r + must be proportional to: the density u, i.e. the more free cadherins are present at a given place, the more of them are susceptible to adhere on the substrate ; the density ρ -v, i.e. the density of free targets available for fixation (this term ensures that v remains bounded by ρ) ; the binding function F (v). We recall that F (v) must be an increasing, positive and bounded function of v. We assume that it depends on the adhesion rate a for cadherins, i.e. the probability of a free cadherin to link to the substrate without any other biological factor. Moreover, the probability a is augmented by the fact that locally other cadherins are fixed, that is the larger is v, the larger should be the value of F (v), and the closer to 0 is v, the smaller should be F (v). We note that, although F (v) is increasing in v, we ask it to be also bounded, so that there is a saturation effect when v gets to its maximum value. Whereas, when v is small, the probability of a free cadherin to bind is small, but it may double when v doubles, and triple when v triples, and so on, so that there F (v) behaves linearly for small v. Since F (v) must be bounded, this linear growth has to saturate for large v. This is a hyperbolic growth behavior and a way to analytically describe it is by means of the tanh(v) function.

Analogously the unbinding rate r -must be proportional to: the density of fixed cadherins v; the unbinding function G(v). We recall that, since local cooperation of bonds strengthen the aggregates, the function G is decreasing, positive and bounded with respect to v. That is the smaller is v, the larger is the value of G(v); and the closer v is to ρ, the smaller is G(v). Note that the density of free cadherins u doesn't play any role in the unbinding rate, so that r -= r -(v). We can then define the unbinding function G(v), analogously as we have done for the binding function F (v), by means of a sigmoidal function.

By the previous considerations, we define F (v) and G(v) as follows,

F (v) = a + tanh(v) a + tanh(ρ) , G(v) = 1 -tanh(αv), (2)
where a > 0 is the adhesion rate and the sigmoid tanh(v) represents the aggregation effect we have previously described. In the definition of the function G(v), the parameter α permits to concentrate or expand the effect of the unbinding : the larger α is, the stiffer the slope of the sigmoid G(v) will be and the narrower will be the effect of the unbinding. A graphic representation of these functions is given in Fig. 2. We also note that the hyperbolic behavior for F (v) and G(v) may be described by other hyperbolic functions, but, with our choice, the behavior (aggregation or not) of the solution to (1) depends only on the two parameters a and α.

Finally, we prove in the following proposition that assuming that the density v is positive and bounded, then both the binding and unbinding functions are bounded, too.

Proposition 1 If v ∈ [0, ρ] then F (v) and G(v) belong to ]0, 1].
Proof. Note that F (v) is a continuous and strictly increasing function on the bounded domain [0, ρ]. Hence its minimum is taken at v = 0, that is F (0) = a/(a + tanh(ρ)) > 0, and its maximum is taken at 2 andα = 1, 2, 3, 4. An increasing parameter a leads to a larger value of the probability function F (v) for small v, that is a larger probability that cadherins bind when (almost) no other links exist; the maximum value for F (v) being always the same this gives a smaller slope in the curve of F (v). When increasing α, the slope of the unbinding probability function becomes larger, so that for large α and v the unbinding is the more and more difficult.

v = ρ, that is F (ρ) = 1. Hence, 0 < F (v) ≤ 1 for all v ∈ [0, ρ]. Concerning G(v), it is a continuous and strictly decreasing function on the bounded domain [0, ρ]. Hence its minimum is taken in v = ρ, that is G(ρ) = 1 -tanh(αρ) > 0, and its maximum is taken in v = 0, that is G(0) = 1. Hence, 0 < G(v) ≤ 1 for all v ∈ [0, ρ].
We finally conclude the description of our model defining the reaction term r(u, v) as follows,

r(u, v) = r + (u, v) -r -(v) = u(ρ -v)F (v) -vG(v), (3)
and endowing the system (1) by suitable smooth initial conditions :

0 ≤ u 0 (x) and 0 ≤ v 0 (x) ≤ ρ (4)
as well as a Neumann homogeneous boundary conditions :

∇u • ν = 0 on ∂Ω (5)
where ν is the exterior unit vector normal to ∂Ω. Note that there is no need to define a boundary condition for v. These boundary conditions ensure that no mass u is lost at the boundary.

Since the global population of cadherins on the membrane is made of free and fixed ones, we can consider that we have a total density distribution of cadherins given by u + v, which mass should be conserved when both u and v evolve. Note that if the initial data (u 0 , v 0 ) is such that:

Ω u 0 (x) + v 0 (x) dx = C, (6)
and since the total mass is conserved, thanks to (5), we have for all time t ≥ 0,

Ω u(x, t) + v(x, t) dx = Ω u 0 (x) + v 0 (x) dx = C.
In order to simplify, in the sequel, we pose C = 1 and the size of the domain |Ω| = 1. Note that the value of the density ρ can vary: if it is larger than 1, then it means that there are much more targets than cadherins so that all cadherins may bind on the substrate in a finite time; if it is smaller than one, then there are not enough targets and only a small proportion of cadherins will bind; if it is equal to 1, then eventually all cadherins may bind, but we will show in simulations that this is not the case.

Analytical results

Now that all the elements of the model are given, in this section, we consider some more theoretical aspects. From a biological point of view, it is clearly important to prove that system (1) admits positive and bounded solutions, otherwise the model is nonsense. Nevertheless, these aspects are beyond the scope of this paper and involve rather technical mathematical theories addressed in a more general framework in [?].

We deal here with the stability of the stationary solutions of the associated problem. Indeed, stable steady states are interesting from a biological point of view, because they characterize the behavior of the studied system at equilibrium, so that the solution of the system should converge for large times to the steady state. As it is known, when considering reaction-diffusion systems in which the ratio of the diffusion coefficients is very small (one of the population diffuses much more slower than the other one), it is possible to obtain stationary solutions which aren't homogeneous in space and which form structures, see [?] for an overview of this kind of models appearing in biology. We note that we are considering an extreme case with respect of those considered in [?], since one of the diffusion coefficient is zero, and we then have a degenerate reaction-diffusion system.

We first consider the solutions to the stationary problem associated to (1), and then pass to the analysis of the steady states of (1), and of their instability, proving the creation of spatial patterns for some choices of the parameter a and α. The stationary problem associated to (1), reads:

∆u = 0 r(u, v) = 0 (7)
completed with Neumann boundary conditions (5) on ∂Ω.

Concerning the existence of a stationary solution to (7), we have the following result.

Proposition 2 Under the hypothesis (6), there exists at least one homogeneous solution (U, V ) to (7) endowed by (5). Moreover, this stationary solution (U, V ) is such that 0 < U < 1 and 0 < V < ρ.

Proof. It is easily seen that the only solution u(x) satisfying both the Neumann boundary condition and the first equation in ( 7) must be constant in x, u(x) = U for all x ∈ Ω. Hence, replacing u(x) = U in (3), it is clear that, if it exists, the function v(x) such that r(U, v(x)) = 0 must be constant too, that is v(x) = V for all x ∈ Ω. Thus, if it exists, the solution (U, V ) to ( 7) is homogeneous in space. Considering now hypothesis (6), with C = 1 and |Ω| = 1, then U = 1 -V . Therefore, the resolution of ( 7) is reduced to find V ∈ [0, ρ] such that:

(ρ -V )(1 -V )F (V ) -V G(V ) = 0.
Let us define the function f : R → R as:

f (v) = (ρ -v)(1 -v) a + tanh(v) a + tanh(ρ) -v(1 -tanh(αv)),
and note that f (0) = aρ/(a+tanh(ρ)) > 0 and that

f (1) = -(1-tanh(α)) < 0.
Moreover, since f is a continuous function on R, then there exists at least one

V ∈]0, 1[ such that f (V ) = 0. Finally, since V ∈]0, 1[, then U = 1 -V ∈]0, 1[.
To conclude the proof we must prove that V is always smaller than ρ.

If ρ > 1, since ]0, 1[⊂]0, ρ[, then V ∈]0, ρ[. Otherwise, if ρ ≤ 1, since f (ρ) = -ρ(1 -tanh(αρ)
) < 0, by the same arguments as before we can conclude that there exists at least one

V ∈]0, ρ[ solving f (v) = 0.
Due to the non linearity of F (v) and G(v), the value V cannot be given analytically, but we can compute it numerically for each suitable choice of parameters a and α. Let us note that, for large values of α, it may become difficult to compute it numerically because there may be two values v both close to 1 and satisfying f (v) = 0. On the other hand, for small values of a, the function f (v) has almost a flat profile when v is close to ρ.

Once proved the existence of the homogeneous stationary solution (U, V ) we can study its stability. We have the following result.

Proposition 3 Under the hypothesis (6), if ln F G V - ρ V (ρ -V ) ≤ 0. (8)
holds, then the homogeneous stationary solution (U,V) is stable. Otherwise, (U, V ) is unstable and patterns formation occurs.

Proof. In order to study the stability of (U, V ) we first have to linearize system (1). Let us note that U and V are linked by the following relation:

U = V G V (ρ -V )F V , (9)
with

F V = F (V ) and G V = G(V ).
Moreover, let us define A = ∂ u r| (U,V ) and B = ∂ v r| (U,V ) . Then a simple computation gives:

A = (ρ -V )F V , B = V (ρ -V )(G V F V -F V G V ) -ρF V G V (ρ -V )F V , (10)
where

F V = dF dv (V ) = 1 -(tanh(V )) 2 a + tanh(ρ)
and

G V = dG dv (V ) = -α(1 -tanh(αV ) 2 ).
Then the linearized system reads:

∂ t u = σ∆u -A u -B v ∂ t v = A u + B v (11)
Let us now seek for a solution with u and v respectively given by:

u = e λt û , v = e λt v, (12)
where û = û(k) and v = v(k) are the respective Fourier transforms. Replacing ( 12) in ( 11) and dividing by e λt , we obtain the following linear system:

λu + σ k 2 u + A u + B v = 0 λv -A u -B v = 0 (13)
Clearly, u = v = 0 is a solution to ( 13), but we seek for non trivial solution to system (13). Hence the determinant of the associated matrix:

M = λ + σ k 2 + A B -A λ -B ,
must be zero. This leads to the following definition of the eigenvalues:

λ = -(σ k 2 + A -B) ± (σ k 2 + A -B) 2 + 4 σ k 2 B 2 , (14)
which sign depends on the sign of A and B.

Note that, in order to have stable solutions, the eigenvalues λ must be negative, so that as time goes to ∞ the solution (u, v) converges to the stationary value (U, V ), see [?] for more details. It is easily seen that A > 0 and that the sign of B depends on the sign of the numerator in (10). If B < 0, then (σ k 2 + A -B) > 0 and 4 σ k 2 B < 0, so that the numerator in ( 14) is always negative, and both eigenvalues are negative too. In this case the stationary solution (U, V ) is stable, and no pattern formation is possible. If B = 0 then the eigenvalues are given by 0 and -2(A+σ k 2 ) which is negative for all k. Thus, in this case, the stationary solution (U, V ) is also stable.

Otherwise, if B > 0 then we have two possibilities. Or 0 < B < A + σ k 2 , or 0 < A + σ k 2 < B. It is easily seen that both cases lead to the same result, that is : one negative and one positive eigenvalue. So that pattern formation may occur.

To conclude the proof, we finally have to show that B ≤ 0 if and only if (8) holds. From (10) we get that B ≤ 0 if and only if

V (ρ -V )(G V F V -F V G V ) ≤ ρF V G V ,
that is, since V , ρ -V , F (v) and G(v) are positive:

G V F V -F V G V F V G V ≤ ρ V (ρ -V )
.

It is now easy to see that the left hand-side in the previous inequality is equal to the first term in (8), concluding the proof.

Due to the complexity of relation (8), it is not possible to give analytical conditions on the parameters a and α such that (8) holds or not. We thus use numerics and show in Fig. 3 that the sign of B is negative for a large range of a and α (the black region). In particular, when α = 1, showing that binding and unbinding forces due to aggregation cannot be symmetric in v if we want patterns formation to occur. On the contrary, the sign of B is positive for a and α in a small domain, the white region. For parameters a and α in this region, the stationary state (U, V ) is thus unstable and patterns formation occurs (see Fig. 4 and6). In order to show the patterns formation in numerical simulations, we thus choose the parameters a and α in the white domain.

Numerical results

In this section we first present some numerical tests showing patterns formation for some theoretical data. Then we take into account the biological aspect of our problem and show on numerical results the different behaviors of the solutions in the biological framework.

We apply a finite differences discretization of system (1). As it is rather classic we just resume briefly our scheme. We define the space step ∆x and ∆y, and the grid points x i = i∆x and x j = j∆x for i, j = 0, 1, 2, . . . N , where N is the number of discretization points. Let t n be the discretized time defined by t n = n∆t, where the time step ∆t must satisfy a stability condition which is given by the discretization of the diffusion term (see ( 15)). Define u n i,j = u(t n , x i , x j ) and v n ij = v(t n , x i , x j ) the approximations of the density functions u and v, and r n ij = r(u n ij , v n ij ) the discretization of the reaction term r(u, v). Then the numerical scheme solving (1) is given by: 0. Determine V and if it does yield to instabilities 2. Initialize v 0 ij and u 0 ij such that (6) holds, then at each time iteration 3. Compute

r n ij = u n ij (ρ -v n ij )F (v n ij ) -v n ij G(v n ij ) 4. Compute u n+1 ij and v n+1 ij by : u n+1 ij = u n ij -σ∆t(Du) n ij -∆t εr n ij v n+1 ij = v n ij + ∆t εr n ij
where (Du) n ij is a classical centered discretization of the second derivatives. 4. Check the stopping criteria Note that the solution u of system (1) will always converge to a homogeneous solution in x, because of the diffusion term, the interesting result is thus the v distribution. Hence, in the following figures we show the density v for a time large enough so that equilibrium is reached.

Another interesting data we can compute is the time evolution of the proportion of fixed cadherins, N v (t), and the one for free cadherins, N u (t). In the sequel we call N v the Immobile Fraction and N u the Mobile Fraction, in agreement with the experimental data showed in Table 1. These two quantities correspond to the zeroth order moments of the solutions v(t) and u(t):

N v (t) = Ω v(t)dx , N u (t) = Ω u(t)dx. ( 16
)
Note also that since we have the total mass conservation, see ( 6), then the evolution of N u (t) is also determined by 1 -N v (t).

These quantities are of interest from a biological and chemical point of view because we know that they must both converge to a constant value (saturation process). In Table 1 we resume some of the corresponding values observed in experiments. Moreover, this convergence is assumed to be exponential and to behave like S -e -kt , where S is the limit saturation value. Then, imposing that N v (t) = S -e -kt , we have that k can be numerically evaluated by means of a linear regression on the function -ln(S -N v (t)). As it will be shown in Fig. 9, the frequency ε plays a central role in order to fit our numerical exponential coefficient k to those obtained in experiments. We finally note that the convergence of N v (t) towards the saturation value S is actually exponential on the first portion of time and slows down when N v gets close to the limit value S.

Patterns formation

Recalling that patterns form when perturbing the homogeneous stationary solution (U, V ), the initial data, in the first two tests, are defined by small perturbations of (U, V ). We define either a random perturbation at each point x ∈ Ω, see Fig. 4, or a perturbation given only in the center x c of the domain Ω, see Fig. 6. These are not biological relevant tests, but we use them to present and discuss the possible behaviors of the solutions with respect to the choice of the parameters a and α. We have fixed the domain Ω to be the square [0, 1] × [0, 1] and we choose N = 100 discretization points on each direction.

In Fig. 4, the perturbation is randomly given on the whole domain Ω. The initial data then reads:

v 0 (x) = V + R 1 (x), (17)
where the perturbation amplitude is represented by R 1 (x) which is a random value of order 10 -3 (left). The patterns formation where particles have aggregated and are fixed to the substrate is clearly shown by the spatial distribution of v (right). Note that at equilibrium the adhesion regions (the set of the lightgray domains) have all the same maximum value max(v) around 0.76, and that the transition to the non-adhesion regions (dark domains) is very sharp. The distribution for the diffusing cadherins u is homogeneous in space, as expected (result not shown).

In Fig. 5 we show the time evolution of N u (t) (dashed line) and N v (t) (continuous line) starting from the initial data (17). After a transition period corresponding to the period during which both functions u and v seems to keep close to the stationary values U and V , a sharp transition takes place at a critical time t c (which value is close to 90 for this run) after which N u and N v converge to a new value. The critical time t c represents the moment Fig. 4 The initial data v 0 given as a random uniform perturbation of (U, V ), see ( 17) (left), and the density v space distribution at equilibrium (right). Note that the two scales in the plots are not the same: v 0 is plotted in the range min(v 0 ) to max(v 0 ), otherwise on 0 to 1 the plot would just be a uniform gray, and v is plotted on the range 0 to 1. Final time T = 500, a 0 = 0.005, α = 1.8, ρ = 1, σ = 0.005 and ε = 1. Fig. 5 The time evolution of Nu(t) (dashed line) and Nv(t) (continuous line) starting from a random perturbation of the stationary solution (U, V ).

at which the structures become spatially delimited. Note that this example is not biologically relevant, since the Immobile Fraction N v doesn't exponentially grow to a saturation level, as it was biologically observed (see [?]). This is due to the fact that the percentage of fixed cadherins at initial time is already higher than the one for free cadherins. Nevertheless, this example numerically confirms that stationary solutions are not stable for a and α opportunely chosen.

Next we consider an initial condition v 0 = v(x, 0) defined by a small gaussian perturbation in the center of the domain Ω, see Fig. 6 (top):

v 0 = V + 0.001 exp -(x -x c ) 2 0.0001 , ( 18
)
with V the stationary value and x c the center of the computational domain, x c = (0.5, 0.5).

Fig. 6 shows the initial data corresponding to ( 18) and the distributions of v when the equilibrium is reached for two different sets of parameters a and α. Note the formation of circular structures where the particles are aggregated.

We did expect the central adhesion region, corresponding to where the gaussian is defined, but other adhesion regions are created around it. This may be caused by the fact that the cadherins aggregation around the gaussian center x c leads to a variation of the densities u and v in the neighborhood of the gaussian, yielding a new perturbation of the stationary value V and thus the formation of new aggregation junctions. This was not the case in the previous test, because the stationary value V was perturbed everywhere in the domain Ω and not enough space was left to other adhesions regions to appear. Again, as expected, a small perturbation of the stationary state, leads to the formation of structures when the parameters a and α are opportunely chosen. The main differences in these two results are the maximum values of v and the size of the circular structures: the larger is α, the larger is the maximum of v (lighter gray) and the larger is the size of the adhesion circular structure.

These numerical tests show that the form and location of the adhesion structures depend on the initial data that we have. In the next section, in order to overcome this problem and to obtain numerical distributions visually more similar to the experimental ones, that is small circular adhesion regions (eventually connected), we construct a more biological relevant initial data and study the influence of the parameters on the equilibrium solution.

Biological framework

We now consider some more biological relevant frameworks and tests. Since our model is defined at a macroscopic scale, it is non-sense to describe the evolution of the densities at the cadherins scale, and thus we assume the domain Ω to be a square of length equal to 10 µm, with N x = N y = 100 grid points, that is a mesh size ∆x = ∆ y = 0.1 µm. We fix the target density value ρ = 1, and, accordingly to [?], we take the diffusion coefficient σ = 3.3 • 10 -2 µm 2 /s. The frequency ε is first fixed to 1 s -1 , and it will be fitted so as to obtain a good agreement for the exponential coefficient k in a second test.

Concerning the initial condition v 0 , we recall that it heavily influences the final distribution v. Since no data is given by experiments, we assume that initially there are no fixed cadherins, so that v(0, x) must be close to zero for all x ∈ Ω. Moreover, we take into account the fact that the membrane may not have a flat profile by defining the initial data v 0 by the sum of N g gaussians randomly distributed on the domain. The choice of each gaussian center is done in such a way that gaussians are sufficiently separated, so that each gaussian may be also seen as a single beginning of aggregate, but these centers are not too far, so that aggregates may also interact. The gaussians eights are also randomly defined while their standard variations are fixed in such a way that each gaussian support covers more than one mesh and may have a small overlap with neighboring gaussians. The initial data v 0 then reads:

v 0 (x) = Ng k=1 R k exp (x -x k ) 2 0.0225 , (19)
where R k is a random value between 0 and 0.005, and the centers x k are such that the distance x k -x i between x k and all others centers x i , with i = k is larger than 0.5. The number of gaussian N g has to be fixed, and we choose it large enough so that it is almost not possible to have N g +1 gaussians satisfying the previous requirements (N g = 266 in our computations). We show in Fig. 7 Fig. 7 The initial data v 0 defined by the sum of random gaussian distributions as in ( 19).

the 2D representation of the initial gaussian distribution we choose. Note that the scale of the plot is between 0 and the maximum value of v 0 , that is 0.005.

Since computational times may be long, we define a stopping criteria by a control on the evolution of the discrete zero moment for the fixed cadherins population, N v . Each niter iterations we compute the Immobile Fraction N v (t niter ) and we compare it to the previous one. When the absolute error between this two quantities is close to zero (e.g., smaller than a given precision 10 -6 ), we stop the computation. Some of the numerical results obtained are given in the following figure.

In Fig. 8 we show the final density v distribution (left) and the time evolution for the Mobile and Immobile Fractions N u and N v (right), for various couples of parameters a and α, and with ρ = 1 and ε = 1. Computations are always stopped when the stopping criteria is satisfied. In A) we choose a small adhesion rate a = 0.001 and a large value of the slope α = 2, the result being a set of interconnected circular regions of adhesion with a maximum value of v around 0.8 and a minimum at almost zero. The Immobile Fraction N v has grown from 0 to almost 0.6 (dot-dashed line), while the Mobile Fraction N u has decreased to almost 0.4 (dashed line). When increasing the adhesion rate to a = 0.05, keeping α = 2, see B), the adhesion regions almost cover the whole domain, yielding to a smaller maximum for v and a larger minimum, but still in a range of values implying that a larger number of cadherins have bound in the adhesion region and a smaller number of free cadherins exists elsewhere. As for test A), we end up with more fixed cadherins than free one, and the values of the Immobile Fraction N v and the Mobile Fraction N u are analogous. Nevertheless, the convergence to the saturation value S is much faster, in fact the Immobile Fraction N v slope in B) is sharper than in A). This is the consequence of the larger adhesion probability a. The same behavior can be observed between C) and D). Note that in these two cases, D) has the smaller value for a, and the main difference with A) and B) stands in the smaller slope value α = 1.6. In both cases the Immobile Fraction N v tends to stabilize around smaller values, 0.5 for C) and 0.4 for B). The final time needed to reach the saturation value S is smaller than for A) and B), while the time needed to reach the equilibrium (stopping criteria) is slightly larger for C) and D), than for A) and B).

When comparing A) with D) and B) with C), that is when fixing the adhesion rates a and changing the slopes α, we note that the maximum value for v is smaller in D) and C) than in A) and B), and also that the minimum is larger. Thus, the influence of the slope α in the unbinding function is translated in the obtention of more interconnected adhesion regions for large α as well as higher values for v. Concerning the evolution of N v and N u , the convergence seems to behave similarly for A) and D), and for B) and C).

We finally note that, although close saturation levels are obtained for the parameters of test A), B) and C), D), the values of the maximum and minimum of the density v may significantly vary. In particular, for both C) and D) the maximum of v is between 0.5 and 0.7, and in C) the minimum is slightly larger than 0 so that we cannot conclude that in the non-adhesion region (dark one) there are no fixed particles. Although we do not have clear biological insight of the real value of the density in the adhesion area, we argue that in the adhesion regions, at equilibrium, the density v may be close to ρ or at least larger than 0.7.

All these examples show a saturation of both the values N v (t) and N u (t) as time is large enough. This fact has been also highlighted in several experiments. We resume in Table 1 the values for the Immobile Fraction which were observed. Note that these values range in the interval 0.3 to 0.7, which also corresponds to what we observe numerically for the N v (t) values. Moreover, the Mobile Fraction in experiments is equal to 1 minus the Immobile Fraction value, and this is the same in our model, since, as already underlined, N u (t) = 1 -N v (t).

We consider now the experiments from [?], for which we know that the saturation value of the Immobile Fraction N v is close to 0.35, and that the convergence has an estimated exponential coefficient K = 13 h -1 = 3.61 • 10 -3 s -1 (see the introduction of Sect. 4 for more details on the exponential behavior). We consider an initial data defined as in ( 19), but such that N v (0) ≈ 0.1, so that the first transitional period, during which the first adhesions appear, reduces to almost zero, but still the Immobile Fraction N v is much smaller than the Mobile Fraction N u . Not considering this choice would just skip the exponential growth of N v (t) for larger times. Moreover, we fix the parameters as follows a = 10 -3 and α = 1.65, ρ = 1, σ = 3.3 • 10 -2 , so that we obtain a saturation level S close to 3.7 and values for the maximum and minimum of the density v respectively close to 1 and 0. The parameter playing a central role here is the frequency ε, varying so as to fit the exponential coefficient k to the experimental one K.

The results of our numerical simulations are shown in Fig. 9, where we plot the N v time evolutions for three values of frequency ε = 0.01, 0.02, 0.03 s -1 . The larger the frequency ε is, the sharper is the slope of the N v curve, that is the smaller is the time needed to reach saturation. Results for some more choices of the relaxation parameter ε, and for the same parameters as in Fig. 9 are given in the Table 2. The coefficient k is computed by means of a linear regression on the curve -ln(S -N v (t)), where S is the saturation limit and time t is restricted to the same initial time interval for all choice of ε. We can conclude that for a fixed set of parameters a, α, ρ and σ the choice of ε significantly influences the exponential rate of convergence, and that in our particular case the best choice is given by ε = 0.02 (note the good agreement of the computed k and the value K obtained in experiments).

We finally consider a last biological situation in which the maximum of the unbinding function G(v) is smaller than 1 so that its efficacy is reduced. It has been shown in [?] that plating cells on antibody coated surfaces did not lead to the formation of aggregates although an even larger fraction of cadherins did bound. Our goal is to numerically prove that we can obtain similar results applying the proposed model. Since we ask for a smaller maximum for the unbinding probability, we slightly modify the definition of the reaction term r(u, v) by multiplying the unbinding function G(v) by a factor γ < 1. Note that this is the same as multiplying by two different efficacy parameters ε the gain and loss terms of the reaction term r(u, v). We shall consider the initial data as defined by ( 19), see also Fig. 7. In Fig. 10 we show the distribution of the density v when the equilibrium is reached in both cases when γ = 1 (left) and γ = 0.5 (right). We note that in the last case, there is only a uniform adhesion region (as it is also mathematically predicted since (8) holds.) Fig. 11 Left: time evolution for the Immobile Fraction Nv(t) for both γ = 1 (continuous line) and γ = 10 -3 (dashed line). For small γ the convergence towards a saturation level is faster, and the saturation level is larger. Right: time evolution of the max(v(t)), the maximum approaches the value of the stationary solution V when γ is small. The parameters are ε = 0.02, ρ = 1, σ = 3.3 • 10 -2 , a = 10 -3 , α = 1.65 and γ = 1 or 0.05, T = 4000.

In Fig. 11, for both γ = 1 (continuous line) and γ = 0.5 (dashed line), we plot on the left the time evolution for the Immobile Fraction N v (t) and on the right the time evolution of the maximum of the density v. We observe that when γ = 0.5 the Immobile Fraction N v has a larger saturation value S than in the case of γ = 1. This is because all the domain Ω is occupied by bound cadherins, and as it can be seen in the right plot the density value for v is almost the same as in the case of γ = 1, as it was observed in experiments, too, see [?].

Conclusion

We propose and study a phenomenological mathematical model intended to describe the cadherins aggregates formation at the cell-cell contacts. This model is composed of two partial differential equations coupled by means of a reaction term which describes the way a cadherin may bind or unbind to targets distributed on a substrate representing a cell membrane, and which takes into account the effects of neighboring pre-existing cadherin-cadherin trans-bonds. Each equation describes the time-space evolution of a population of cadherins: the free ones, not linked to the substrate and diffusing, and the fixed ones, not diffusing. After having numerically proved that the model is able to reproduce the formation of aggregates, we construct a more biologically relevant framework and compare the results obtained by means of numerical simulations to experimental ones. In particular, the proposed model shows a saturation of the density of both populations which is of the same order of the one observed in several experiments, see Table 1 and references therein. Moreover, suitably choosing the parameters of the model we are able to compute the supposed undergoing exponential convergence coefficient with a very good agreement with the one observed in [?]. Finally, adapting the efficacy parameter ε in the gain and loss parts of the reaction term r(u, v), in such a way that the loss term becomes negligible, the model justifies the uniform adhesion versus aggregates adhesions as experimentally shown in [?]. This work is a first attempt to model the cell-cell junctions creation at a macroscopic scale, and may be enriched and specified including more biophysical aspects, such as cadherin cis-interactions, or cadherin association to F-actin including the retrograde effect of the actin filaments. Moreover, this work may have a much wider implication in biological processes since many membrane-born proteins associate in clusters at the cell surface. Their recruitment in discrete aggregates is often associated to their specific function and understanding the biochemical and physical principles that govern this process is of major importance.

Figure 1 :

 1 Figure 1: (A): example of accumulation of immunofluorescently labelled cadherin adhesion complexes in cells seeded on a glass coverslip coated at saturation with recombinant cadherin ectodomain. The signal is detected by total internal reflection fluorescence microscopy (TIRFM); the brightest is the signal, the stronger is the local accumulation of cadherins. Notice the accumulation is hot spots. (B) : Reflection interference contrast microscopy (RICM) imaging of the same cells objectivizing the distance (h) separating locally the cell membrane from the glass coverslip. The closest the membrane is from the membrane, the darker is the signal. Scale bar: 5 µm. Notice the preferential accumulation of cadherins in areas of the lowest h value. The aligned hot spots at the cell periphery are associated to radial actin cables pushing the membrane close to the substratum. Extracted from [Lam+07].

Fig. 1

 1 Fig. 1 (A): example of accumulation of immunofluorescently labelled cadherin in cells seeded on a glass coverslip coated at saturation with recombinant cadherin ectodomains. The signal is detected by total internal reflection fluorescence microscopy (TIRFM); the brightest is the signal, the stronger is the local accumulation of cadherins. Notice the accumulation in hot spots. (B) : Reflection interference contrast microscopy (RICM) imaging of the same cell region underlying the distance separating locally the cell membrane from the glass coverslip. The closest the membrane is from the membrane, the darker is the signal. Scale bar: 5 µm. Notice the preferential accumulation of cadherins in areas of lowest membrane distance. The aligned hot spots at the cell periphery are associated to radial actin cables pushing the membrane close to the substratum. Extracted from [?].

Fig. 2

 2 Fig. 2 Examples of the probability functions F (v) (left) and G(v) (right) for different values of the parameters : a = 0.025, 0.05, 0.1, 0.2 and α = 1, 2, 3, 4.An increasing parameter a leads to a larger value of the probability function F (v) for small v, that is a larger probability that cadherins bind when (almost) no other links exist; the maximum value for F (v) being always the same this gives a smaller slope in the curve of F (v). When increasing α, the slope of the unbinding probability function becomes larger, so that for large α and v the unbinding is the more and more difficult.

Fig. 3

 3 Fig.3The sign of B in the parameters domain (a, α) =]0, 1] × [0, 3]. Black region corresponds to B ≤ 0 (stable solution (U, V )). White region corresponds to couples of parameters a, α for which B > 0 (unstable solution (U, V )).

Fig. 6

 6 Fig. 6 Top: initial data v 0 given by a small gaussian perturbation of the stationary value V : small white point in the middle of the domain. Bottom: density distributions for v at equilibrium (final time T = 1500). Left: the solution for the parameters a = 0.005, α = 1.8, ρ = 1 and σ = 0.05. Right: the solution for the parameters a = 0.005, α = 1.4, ρ = 1, σ = 0.05 and ε = 1.

Fig. 8

 8 Fig.8The final 2D representation density v (left) and the time evolution for Nu and Nv. The parameters are ρ = 1, σ = 3.3 • 10 -2 , ε = 1 and for figures: A) a = 0.001, α = 2; B) a = 0.005, α = 2; C) a = 0.005, α = 1.6; D) a = 0.001, α = 1.6.

Fig. 9

 9 Fig. 9 Time evolutions of Nv for different values of the frequency : ε = 0.01 dot-dashed line, ε = 0.02 continuous line, ε = 0.03 dashed line. The parameters are ρ = 1, σ = 3.3•10 -2 , a = 10 -3 , α = 1.65.

Fig. 10

 10 Fig. 10 Left: final equilibrium v distribution (T = 4000) for γ = 1, we note the formation of adhesion regions. Right: final equilibrium v distribution (T = 4000) for γ = 0.5, the distribution is homogeneous in space. The parameters are ε = 0.02, ρ = 1, σ = 3.3 • 10 -2 , a = 10 -3 , α = 1.65 and γ = 1 or 0.5.

Table 1

 1 Experimental saturation values for the Immobile Fraction, corresponding to N (v) at equilibrium.

	Cadherins type Cell Type	Immobile Fraction Mobile Fraction Reference
	E-cadherin	MDCK	0.6-0.7	0.3-0.4	[?]
	N-cadherin	C2C12	0.6	0.4	[?]
	N-cadherin	neurons	0.3	0.7	[?]
	E-cadherin	A431	0.3-0.4	0.6-0.7	[?]
	E-cadherin	MDCK	0.5-0.6	0.5-0.6	unpublished
	N-cadherin	HEK	0.4-0.6	0.4-0.6	unpublished

Table 2

 2 Relaxation parameters τ and corresponding exponential coefficients

Acknowledgements POS has been supported by an ANR funding (2010 Blan 1515). BS is supported by the Labex LMH (ANR-11-IDEX-003-02) and by the FMJH (ANR-Investissements d'Avenir). The authors thank the CNRS, INRIA and INSERM for the founding of the PEPS-MBI project "MAC: Modélisation d'Adhésion des Cadhérines" the ANR 2010 Blan 1515 and Human Frontier Science Program grant RTG0040/2012, and the ANR blanche project "Kibord" (ANR-13-BS01-0004) funded by the French Ministry of Research.