Erin R Higgins

Heiko Schmidle

Martin Falcke
email: falcke@hmi.de

Waiting time distributions for clusters of IP3 receptors

Keywords: Inositol 1, 4, 5-trisphosphate receptor, calcium dynamics, ion channels 1

HAL is

A c c e p t e d m a n u s c r i p t 1 Introduction

The inositol 1, 4, 5-trisphosphate receptor (IP 3 R) acts as a channel which controls release of Ca 2+ from the endoplasmic reticulum (ER) into the cytosol of a cell. The al. [START_REF] Rüdiger | Hybrid Stochastic and Deterministic Simulations of Calcium Blips[END_REF][START_REF] Shuai | A Kinetic Model of Single Clustered IP 3 Receptors in the Absence of Ca 2+ Feedback[END_REF] (model 1), the same model modified to include sequential binding of IP 3

and Ca 2+ as suggested by Adkins and Taylor [START_REF] Adkins | Lateral Inhibition of inositol 1,4,5trisphosphate receptors by cytosolic Ca 2+[END_REF] (model 2), and the model developed by Sneyd and Dufour [START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF] (model 3) as given in reference Ullah et al. [START_REF] Ullah | Modeling the Statistics of Elementary Calcium Release Events[END_REF]. The IP 3 R is comprised of four subunits, which are assumed to be identical, and these three models each describe the transitions between subunit states, most of which result from the binding and release of IP 3 R and Ca 2+ . Model 1 [START_REF] Shuai | A Kinetic Model of Single Clustered IP 3 Receptors in the Absence of Ca 2+ Feedback[END_REF] is closely based on the DeYoung-Keizer model [START_REF] De Young | A Single-pool inositol 1,4,5-trisphosphatereceptor-based model for agonist-stimulated oscillations in Ca 2+ concentration[END_REF], and assumes that a subunit has binding sites for IP 3 , activating Ca 2+ and inhibitory Ca 2+ . The model by Sneyd and Dufour [START_REF] Ullah | Modeling the Statistics of Elementary Calcium Release Events[END_REF][START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF][START_REF] Falcke | Reading the patterns in living cells -the physics of Ca 2+ signaling[END_REF] was designed to agree with several experimental observations.

While it has been known for about a decade now that puffs are random release events, we recently showed that spike sequences of intracellular Ca 2+ oscillations are random events also [START_REF] Skupin | How does intracellular Ca 2+ oscillate: by chance or by the clock?[END_REF]. The random opening and closing of individual channels, i.e.

A c c e p t e d m a n u s c r i p t

the random state transitions of their subunits, is believed to be the most important source of noise causing the random behavior of the whole system [START_REF] Falcke | On the role of stochastic channel behavior in intracellular Ca 2+ dynamics[END_REF]. That requires stochastic theory across all structural levels from a channel subunit up to cell level for modelling. However, existing subunit models typically comprise 6 or more subunit states, leading to hundreds of channel states and many thousands of channel cluster states. That makes stochastic theory on cluster level already almost impossible, not to mention cell level. Thus theoretical studies on puff characteristics used either extremely simplified channel models (with two subunit states [START_REF] Shuai | Stochastic properties of Ca 2+ release of inositol 1,4,5trisphosphate receptor clusters[END_REF], two or four channel states [START_REF] Nguyen | A stochastic automata network descriptor for markov chain models of instantaneously-coupled intracellular Ca 2+ channels[END_REF]) or relied on simulations [START_REF] Ullah | Modeling the Statistics of Elementary Calcium Release Events[END_REF].

Recently we presented a method which circumvents writing down the master equation of the channel molecule with its huge number of states and gets by with the subunit master equation only in calculations of waiting time distributions for channel clusters [START_REF] Thul | Waiting time distributions for clusters of complex molecules[END_REF]. Here we follow a different approach. Channel states are defined by the number of subunits in subunit states. We derive the channel Master equation from the subunit Master equation. Subsequently, we lump channel states into aggregates based on time scale separation. Standard methods for calculating cluster properties can then be applied to the Master equation defined for the lumped states.

Using these approximations, we compared the probability distributions for the first opening time of a cluster and the first closing time of a cluster, at various IP 3 R and Ca 2+ concentrations. These distributions were calculated analytically as well as through simulations, and assume a cluster consists of five channels. The three models

Single Channel Activation

We would like to calculate F (A, t|I, 0), the probability distribution for a channel first being activated at time t, given it is in state I at time t = 0. A denotes the set of active channel states, and we say a channel is activated when it is in one of these states. We are assuming we have the subunit transition matrix W , where w i,j dt gives the probability for a subunit to move from state x j to state x i in infinitesimal time dt, for i = j. w j,j = -n i=1,i =j w i,j where n is the number of subunit states. The probability for a subunit to stay in state x i during time dt is 1 + w i,i dt.

We demonstrate how to build the subunit transition matrix W with model 1 (see section 5 for a detailed description of the model). First, we order the subunit states:

(110 + Active, 111, 010, 011, 100, 101, 000, 001). For instance the transition rate from state 100 to 110+Active is a 5 c where c denotes the Ca 2+ concentration. Now we can build the 8 × 8 subunit transition matrix

W = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 b 2 a 1 • p 0 a 5 • c 0 0 0 b 0 a 2 a 0 +b 0 • c -2 0 a 3 • p 0 a 5 • c 0 0 b 0 b 1 a 0 +b 0 0 -3 b 4 0 0 a 5 • c 0 0 b 3 a 4 • c -4 0 0 0 a 5 • c b 0 b 5 a 0 +b 0 0 0 0 -5 b 2 a 1 • p 0 0 b 5 0 0 a 2 • c -6 0 a 3 • p 0 0 b 5 0 b 1 0 -7 b 4 0 0 0 b 5 0 b 3 a 4 • c -8 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
All entries have dimension s -1 , p denotes the IP 3 concentration.i represents the negative sum of the ith column.

We can use W to form the channel state transition matrix P . An entry in P , p i,j dt,

gives the probability for a channel to move from state X j to state X i in infinitesimal

A c c e p t e d m a n u s c r i p t

time dt, for i = j. The diagonal elements are given by p j,j = -r i=1,i =j p i,j . The off-diagonal elements of P are found as follows. A channel state is represented by a vector V where an element of the vector, v i , gives the number of subunits in subunit state x i . Thus the length of V is n, and the sum of the elements in V is the number of subunits in each channel, h. The transition probability p i,j is 0, if X j and X i differ by more than one subunit and p i,j = φ w k,l , if X i can be reached from X j by one subunit in state x l moving to state x k , and X j has φ subunits in state x l .

To find F (A, t|I, 0) we first remove all rows and columns in P which correspond to active channel states and denote the new matrix by P . That is, if X i ∈ A then we remove row i and column i. This means if the system reaches an active state, it cannot leave it. Let y(X i , t|I, 0) be the probability of being in nonactive state X i at time t and let y = {y(X i , t|I, 0), X i / ∈ A}. Then y can be found by solving the master equation ẏ = P y.

Let f (A, t|I, 0) be the probability that an activated state is reached in the time interval [0, t], let r be the total number of channel states and let r a be the number of open channel states. Then

f (A, t|I, 0) = 1 - i:X i / ∈A y(X i , t|I, 0) = 1 - r-ra i=1 r-ra k=1 c k V k,i e λ k t
where λ k and V k , k = 1...rr a are the eigenvalues and eigenvectors respectively, of P and where c k , k = 1...rr a are determined by the initial conditions. Then

F (A, t|I, 0) = df (A, t|I, 0) dt = - r-ra i=1 r-ra k=1 c k λ k V k,i e λ k t .

Cluster Activation

A cluster is first activated when any of the channels in the cluster are first activated.

We assume the initial state of a channel is unknown, so we use the weighted average over all possible initial channel states, where the weights are given by the probabilities

of
F open (A, t) = NF ch (A, t)G(A, t) N -1
where

G(A, t) = 1 - t 0 F ch (A, τ )dτ
gives the probability that a channel has not been activated by time t and N is the number of channels in the cluster. That is, the probability for the cluster to be first

activated
F open (A, 0) = 1 -(1 -F ch (A, 0)) N .
That is, the probability that the cluster opens instantly is the negation of no channel opening instantly.

We use the non-open states weighted according to their resting probabilities as initial condition for calculating F (A, t) for t > 0. The resting probabilities for the initial channel states are given by the normalized eigenvector of P corresponding to the eigenvalue λ = 0.

Cluster Closure

The cluster of channels is considered closed when all channels in the cluster are closed. Let F close (A, t) denote the probability distribution for a cluster of channels

A c c e p t e d m a n u s c r i p t

to first close at time t if it has opened at time t=0 (closing time distribution). To calculate F close (A, t), we find the cluster state transition matrix, M, from the channel state transition matrix P , in an analogous way to how we found P from the subunit state transition matrix, W . That is, a cluster state is represented by a vector with length given by the number of channel states. An entry in the vector gives the number of channels in the corresponding channel state. Thus, the sum of the elements in the vector is the number of channels in the cluster, N. The probability of a transition between cluster states occurring which involves more than one channel changing state instantaneously is infinitely small compared to the probability of a transition where only one channel changes state, and so for such transitions, the entry in M is 0. If cluster state X i can be reached from state X j by one channel in state X l moving to state X k , and X j has φ channels in state X l , then the entry in M, m i,j = φ p k,l . The diagonal entries of M are such that the columns of M sum to 0.

Once M has been constructed, F close (A, t) can be found analogously to the way in which F (A, t|I, 0) is found in section 2. That is, let M be the cluster state transition matrix with rows and columns corresponding to closed cluster states removed. Let ρ be the total number of cluster states and let ρ a be the number of closed cluster states.

Let μ k and U k , k = 1...ρρ a be the eigenvalues and eigenvectors respectively, of M .

Then

F close (A, t) = - ρ-ρa i=1 ρ-ρa k=1 d k μ k U k,i e μ k t .
where d k , k = 1...ρρ a are determined by the initial conditions.

We choose as the initial probabilities for a channel to be in each of the channel states the probabilities at the expected cluster opening time. That is, we begin with the resting probabilities, then solve the master equation,

ż(t) = P 0 z(t),
where z(t) gives the probabilities of being in the channel states at time t, and P 0 The large number of cluster states presents a problem when calculating F close (A, t).

The number of channel states, r is given by r = h+n-1 n-1 . Recall, n is the number of subunit states and h is the number of subunits per channel. The number of cluster states, ρ, is given by ρ

= N +r-1 r-1
. Recall, N is the number of channels in the cluster.

So ρ is very large for the values of n, h and N we require, which results in the above method for finding F close (A, t) being too slow and requiring too much memory to be feasible.

To overcome the problem of there being too many cluster states, we reduce the number of channel states before finding M. We developed a reduction method motivated by and similar to that described in Huisinga et al. [START_REF] Deuflhard | Identification of almost invariant aggregates in reversible nearly uncoupled Markov chains[END_REF]. It is based on the theory of weakly coupled Markov chains. The weakly coupled state aggregates are identified by a Perron cluster analysis in the study by Huisinga et al. We replace cluster identification by a simpler method. In one case we compared our results with those of the method used in [START_REF] Deuflhard | Identification of almost invariant aggregates in reversible nearly uncoupled Markov chains[END_REF] and found them to be identical.

We identify aggregates by repeatedly finding the largest transition rate in the channel state transition matrix P , then combining the two states connected by the transition, into one aggregate state. In general we combine two states X i , X j to an aggregate state {X i , X j } and after repeating the reduction process we get a set of aggregate states of the form {X a i , i = 1 . . . m}. The transition rate for moving from the aggregate state

{X a i , i = 1 . . . m 1 } to the aggregate state {X b i , i = 1 . . . m 2 } is then given by m 2 j=1 m 1 i=1 p b j ,a i V 0 (a i ) / m 1 i=1 V 0 (a i ) .
V 0 represents the stationary solution of the whole un-reduced system and p b j ,a i denotes the transition probability from state a i to b j . This is repeated until the number of aggregate channel states is sufficiently small that the number of cluster states is small

Subunit Models

We have compared results using three models of IP 3 channel subunits. These include the nine state model of Rüdiger et al. [START_REF] Rüdiger | Hybrid Stochastic and Deterministic Simulations of Calcium Blips[END_REF] which we have reduced to eight states, the same model but with conditions of sequential binding of Ca 2+ and IP 3 imposed, and the model of Sneyd and Dufour [START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF] modified and reduced to six states as given by Ullah and Jung [START_REF] Ullah | Modeling the Statistics of Elementary Calcium Release Events[END_REF].

Model 1

Model 1 is based on the DeYoung-Keizer model [START_REF] De Young | A Single-pool inositol 1,4,5-trisphosphatereceptor-based model for agonist-stimulated oscillations in Ca 2+ concentration[END_REF], whereby the subunit has an IP 3 binding site, an activating Ca 2+ binding site and an inhibitory Ca 2+ binding site.

Subunit states are represented by the triplet (ijk), where i, j and k represent these three binding sites respectively. An occupied site is represented with a 1 and an unoccupied site with a 0. In the DeYoung-Keizer model [START_REF] De Young | A Single-pool inositol 1,4,5-trisphosphatereceptor-based model for agonist-stimulated oscillations in Ca 2+ concentration[END_REF], the subunit is considered active when it is in state (110). That is, when the IP 3 and activating Ca 2+ binding sites are occupied but the inhibitory Ca 2+ binding site is not occupied. Rüdiger et al. [START_REF] Rüdiger | Hybrid Stochastic and Deterministic Simulations of Calcium Blips[END_REF] have added an extra state labelled 'Active' with transitions to and from state (110). The subunit is only active when it is in this state. This model is illustrated in Fig. 1. We have combined the 'Active' state and state (110) into one state by assuming the transition rates between the two states are fast, and so reduced the model to eight states. We have used the transition rates given by Shuai et al. [START_REF] Shuai | A Kinetic Model of Single Clustered IP 3 Receptors in the Absence of Ca 2+ Feedback[END_REF],

A c c e p t e d m a n u s c r i p t

which are listed in the Appendix.

A channel contains four such subunits, and is considered open when any three of the four subunits are in the active state.

Model 2

It had been suggested by Adkins and Taylor [START_REF] Adkins | Lateral Inhibition of inositol 1,4,5trisphosphate receptors by cytosolic Ca 2+[END_REF] that the binding of Ca 2+ to the activating and inhibitory binding sites may be affected by whether or not the IP binding rate when IP 3 is bound, are made small (they are set to 1.0 × 10 -3 μM -1 s -1 ).

In order to maintain the condition of detailed balance, the corresponding reverse transition rates (the b i values) have also been made small, so that the dissociation constants, K i = b i /a i remain the same. The transitions that are made slow are indicated by dotted lines in Fig. 1.

Model 3

Sneyd and Dufour [START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF] describe a 10 state model of the IP 3 R subunit. It includes one IP 3 binding site, two binding sites for Ca 2+ activation and one binding site for Ca 2+

inactivation. The Sneyd and Dufour [START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF] model was designed to produce agreement with a number of results from experimental data. Falcke [START_REF] Falcke | Reading the patterns in living cells -the physics of Ca 2+ signaling[END_REF] modified the model to overcome a problem involving a lack of Ca 2+ conservation in triangular motifs employed to achieve saturating binding kinetics. These triangular motifs are replaced by square motifs Fig. 2 B, and the full model contains then thirteen states. The resulting transition diagram is shown in Fig. 2. We have used the parameter values given by Ullah and Jung [START_REF] Ullah | Modeling the Statistics of Elementary Calcium Release Events[END_REF], but changed them slightly in order to obey detailed balance. They can be found in the Appendix.

A channel consists of four of these subunits, and is considered open when all four

Results

We show results for a cluster consisting of five channels, using various values for the In Fig. 3 A,B ,we see that the opening distributions for models 1 and 2 are very similar when c close < 0.1 μM. At large Ca 2+ concentrations, model 2 has a much longer expected opening time compared to model 1 (see Fig. 4). At high Ca 2+ concentrations, there is a high probability of a subunit being in state (011) at rest. In model 1, the fastest route from (011) to (110) is through the binding of IP 3 followed by the release of inhibitory Ca 2+ . However, in model 2 this route is not fast, resulting in a longer expected opening time. When the Ca 2+ concentration is low (∼ 0.1 μM) subunits are more likely to initially be in state (000), from which transitions to the active state are very fast in both models. However, if the Ca 2+ concentration is very low (0.01 μM for example), the binding of Ca 2+ to the activating Ca 2+ site is slow, resulting in the increase in the expected opening times seen in Fig. 4 A,B.

In model 3, the single Ca 2+ dependent transition rate into an active state depends on the concentration of IP 3 too, which is slow for small IP 3 concentration. Therefore, the effect of increasing the Ca 2+ concentration is to increase the transition rates into show that the agreement of the analytical and numerical distributions for small times is rather good. However, if the distribution is more spread out, a small error in the tail of the analytic distributions may be causing the bad agreement. . This is the maximum number for which results could be computed due to limitations set by computational efficiency. 1 We see also from Fig. 14 that this is close to the minimum value for which we can expect reasonably accurate results, when using a single channel.

A c c e p t e d m a n u s c r i p t 7 Stop criterion for Reduction

It can be seen in Fig. 14 A are reasonably accurate when using a single channel, they will also be reasonably accurate when using 5 channels. Thus, the method we suggest for predicting the accuracy of results after aggregation, is to compute the expected closing time for a single channel using the full model and to compare this with the expected closing time computed after the number of states has been reduced. The minimum number of aggregate states for which reasonable approximations can be expected is the number just above the dramatic change of results.

A c c e p t e d m a n u s c r i p t 8 Discussion

Many models describing the kinetics of IP 3 R subunits have been proposed. We have compared the probability distributions for the expected opening and closing times of a cluster of five IP 3 Rs resulting from three different subunit models. The first model is that of Shuai et al. [START_REF] Shuai | A Kinetic Model of Single Clustered IP 3 Receptors in the Absence of Ca 2+ Feedback[END_REF], which is similar to the DeYoung-Keizer model [START_REF] De Young | A Single-pool inositol 1,4,5-trisphosphatereceptor-based model for agonist-stimulated oscillations in Ca 2+ concentration[END_REF]. Here there are binding sites for IP 3 , activating Ca 2+ and inhibitory Ca 2+ , resulting in eight subunit states. Shuai et al. [START_REF] Shuai | A Kinetic Model of Single Clustered IP 3 Receptors in the Absence of Ca 2+ Feedback[END_REF] add an extra 'Active' state reachable only from the state where IP 3 and activating Ca 2+ are bound, and inhibitory Ca 2+ is not bound.

We have reduced these two states to one state by assuming the transitions between them are fast. The parameter values are such that the binding of IP 3 reduces the affinity of the inhibitory Ca 2+ site, and Ca 2+ binding to the inhibitory site reduces the affinity of the IP 3 binding site. The second model used is a modification of the first suggested by Adkins and Taylor [START_REF] Adkins | Lateral Inhibition of inositol 1,4,5trisphosphate receptors by cytosolic Ca 2+[END_REF]. Here the inhibitory Ca 2+ binding rate when IP 3 is bound is made very small, as is the activating Ca 2+ binding rate when IP 3 is not bound. The reverse transition rates are also made small to maintain detailed balance. The third model is a thirteen state model given by Sneyd and Dufour [START_REF] Sneyd | A dynamic model of the type-2 inositol trisphosphate receptor[END_REF] and extended by Falcke [START_REF] Falcke | Reading the patterns in living cells -the physics of Ca 2+ signaling[END_REF] designed to accommodate various experimental findings.

We have compared the opening and closing time distributions resulting from the three models, at various Ca 2+ and IP 3 concentrations, as well as the expected opening and closing times. We found the three models gave significantly different results. At low Ca 2+ concentrations (< 0.1 μM), models 1 and 2 have similar expected opening times, but at larger Ca 2+ concentrations model 2 has a much longer expected opening time. When the Ca 2+ concentration is large, both the activating and inhibiting Ca 2+

binding sites are likely to be occupied prior to the addition of IP 3 . Model 1 then opens comparatively quickly due to the binding of IP 3 followed by the release of inhibitory Ca 2+ . This route to activation is slow in model 2 as the rate for the release of inhibitory Ca 2+ when IP 3 is bound is small.

The models show also differences in closing behavior. Using model 1 there is a decrease in the expected closing time of a cluster, as the Ca 2+ concentration increases We have used a fixed value for the Ca 2+ concentration when the cluster is open.

However, the concentration around the cluster may vary depending on whether indi-

vidual channels are open or closed. In Fig. 11 we use a variable Ca 2+ concentration, so that channels in the cluster which are closed see a lower Ca 2+ concentration than

those that are open. Making the Ca 2+ concentration variable has had very little impact. Fraiman et al. [START_REF] Fraiman | Analysis of Puff Dynamics in Oocytes: Interdependence of Puff Amplitude and Interpuff Interval[END_REF] note that the time taken for a single inactivated IP 3 R to become uninhibited is much longer than a typical puff duration. Our results support this statement. Consequently, the Ca 2+ concentration seen by a closed IP 3 R has little effect, as all channels are likely to close before it has time to reopen. These results imply that puffs with these models are insensitive towards spatial arrangement of individual channels within the cluster as long as the distance is not too large. We mimicked a decay of the Ca 2+ concentration on the channel distance corresponding to 3-4.6 diffusion lengths of free Ca 2+ .

Our method for calculating the closing time distribution of a cluster requires the number of channels in the cluster to be small (≤ 5) so that the number of possible cluster states is not too large. There have been various estimates of the number of channels in a cluster. Some estimate it to be in the range of 20-30 [START_REF] Callamaras | Activation and co-ordination of InsP 3 -mediated elementary Ca 2+ events during global Ca 2+ signals in Xenopus oocytes[END_REF][START_REF] Dupont | Quantal release, incremental detection, and long-period Ca 2+ oscillations in a model based on regulatory Ca 2+ -binding sites along the permeation pathway[END_REF][START_REF] Sun | A continuum of InsP 3 -mediated elementary Ca 2+ signalling events in Xenopus oocytes[END_REF], while In order to calculate the closing time distribution of a cluster, we reduced the number of channel states so that the number of cluster states was not too large.

We tried different methods to reduce the system among them a rapid-equilibrium approximation: We identified the states connected by the fastest transition rates, assumed equilibrium just between these two states, lumped them accordingly and

A c c e p t e d m a n u s c r i p t

repeated the procedure with the new set of states. With respect to accuracy, that method was as good as the one we used here. Its only disadvantage was that it did not allow for re-separating closed and open states once the aggregates had been found.

The state aggregation introduces error to the calculation. This is particularly apparent in the poor agreement between the analytic and simulation estimates of the expected cluster closing times for model 2. There is a loose relationship between the expected opening time, and the agreement between simulation and analytic results, such that if the expected opening time is long it is likely the agreement will be worse than if the expected opening time is short. Thus, it seems that the error results from the smallest eigenvalues of the coefficient matrix.

We tried a variety of criteria for the degree of state aggregation still providing good approximations. We required for instance the inverse of the smallest transition rate lumped to be much smaller than moments of the waiting time distributions or residence times of the system in aggregates. These are criteria suggested by the use of time scale separation as the basis for state aggregation. Surprisingly, meeting conditions of that type did not give conclusive predictions on the quality of the approximation. The criterium to stop lumping states into aggregates when analytic single channel results get bad compared to single channel simulations was the most predictive and feasible criterium.

We have not related our results to experimental observations in this study yet.

Fitting models to experimental data will be the subject of future research. Puffs can be found in a very small range of IP 3 concentration only below which no release is observed and above which spiking or sustained release is found. Only recently. the lab of I. Parker's developed a method to observe puffs over a larger range of IP 3 concentrations by decoupling puff sites by the buffer EGTA [22]. This method will produce data with which we can meaningfully compare our results. C: F open (A, t) for model 3.

k 1 2.0 μM -1 s -1 k -1 0.04 s -1 k 2 37.4 μM -1 s -1 k -2 1.4 s -1 k 3 0.11 s -1 k -3 29.8 s -1 k 4 4.0 μM -1 s -1 k -4 0.37 s -1 k 5 2.0 μM -1 s -1 l 1 10.0 μM -1 s -1 l 3 100.0 μM -1 s -1 l 5 0.1 μM -1 s -1 L 1 0.

IP 3

 3 Rs release Ca 2+ in response to an increase in the concentration of inositol 1, 4, 5trisphosphate (IP 3 ) in the cytosol, which results from a stimulus, such as the binding of a hormone to the cell membrane. IP 3 Rs open and close randomly, and the open and closing probabilities are determined by the binding of IP 3 and Ca 2+ . A small increase in the concentration of cytosolic Ca 2+ stimulates channel activation, whereas a large increase causes inhibition of the IP 3 R. This results in a bell-shaped dependence of the open probability on the cytosolic Ca 2+ concentration. Activation of the channels takes place on a faster time scale than inhibition. IP 3 Rs are found in clusters, and so the release of Ca 2+ resulting from the opening of one channel in a cluster, stimulates the opening of other channels in the cluster. The release of Ca 2+ from the cluster is referred to as a Ca 2+ puff. Puffs are the elemental events of intracellular Ca 2+ release. Various models have been proposed for describing the IP 3 R [15, 7]. Here we have compared the opening and closing probability distributions resulting from three different models. These include the model proposed in Rüdiger et al. and Shuai et

A c c e p t e d m a n u s c r i p t

  exhibit significantly different behaviour. The sequential binding of IP 3 and Ca 2+ imposed in model 2 results in longer expected opening times compared to model 1 and longer expected closing times compared to both models 1 and 3. Model 3 has a shorter expected closing time than models 1 and 2, but has a much longer expected opening time at low IP 3 concentrations and high Ca 2+ concentrations. These probability distributions are calculated assuming the Ca 2+ concentration in the vicinity of the cluster remains constant. That is, we use one fixed value for the Ca 2+ concentration when the cluster is closed, and another fixed value for Ca 2+ when the cluster is open, i.e., when at least one channel is open. We have also looked at the effect of using a variable Ca 2+ concentration while the cluster is open, whereby a closed channel sees a lower concentration than an open channel. The resulting probability distributions are very similar to when a fixed value is used, in all three models.

A c c e p t e d m a n u s c r i p t F

 t is found from P by setting transition rates out of open states to 0, because we ave interested in those trajectories only, which did not go through open states. z(t) is then evaluated at the expected opening time for the cluster, which is calculated from open (A, t). Initially, we assume one channel in the cluster is open, while the others are closed. To enforce this condition, we find from z(t) the conditional probability of being in each state where exactly one channel is open, given that the cluster is in one such state. All other cluster states have an initial probability of zero.

A c c e p t e d m a n u s c r i p t

  enough for computations to be feasible or till the stop criterion has been reached (see section 7). This process of combining channel states results in open and closed channel states being combined. However, as we require open states to be distinct from closed states, we divide the aggregate states which contain both open and closed states into two separate states, one containing only the closed states and one containing only the open states. After these aggregate states are split, transition rates to and from the new aggregate states must be found.

A c c e p t e d m a n u

  s c r i p t are in state O , Õ, Ō (the open states) or state Ã, Ā (the activated states) or an intermediate combination, such as one subunit in state O and three in state Ā.

Ca

  2+ concentration when the cluster is closed, c close , the Ca 2+ concentration when the cluster is open, c open , and the IP 3 concentration, p. When finding the steady channel state probabilities, we use p = 10 -7 and the Ca 2+ concentration c close . When the cluster is open, we use the Ca 2+ concentration c open and it remains fixed regardless of the number of open channels. When c close , c open and p are not being varied, their values are 0.1 μM, 100 μM and 0.15 μM respectively. Figures 3-9 show the effect of changing c close , c open and p when using models 1, 2 and 3. We present the probability distribution for the opening time F open (A, t) and closing time F close (A, t) of a cluster of channels and the expected value of the opening time distribution T open and closing time distribution T close .

A c c e p t e d m a n u s c r i p t

  the inactivated or shut states. That slows down activation (Fig.3C, Fig.4) and the expected opening time does not exhibit a minimum in dependence on the Ca 2+ concentration.In Fig.5A-C and Fig.6A,B we see that increasing p results in an initial rapid decrease in the expected opening time. All three models exhibit high sensitivity to the IP 3 concentration at low concentrations, and low sensitivity at high concentrations.Model 2 again has a long expected opening time compared to model 1 as the fastest routes to the active state in model 1 are made slow in model 2.

Fig. 7 AA c c e p t e d m a n u s c r i p t 100

 7100 Fig. 7 A,B and Fig. 8 show that changing the Ca 2+ concentration has a strong effect on the closing behavior of model 1, but has very little effect on the closing time distribution of model 2. In model 2, the Ca 2+ dependent transition rate out of the open state is very low, so the Ca 2+ concentration has little effect on the expected closing time. The most prominent closing transition is IP 3 dissociation. At high Ca 2+ concentrations, model 1 closes quickly due to the binding of inhibitory Ca 2+ . However, in model 2 inhibitory Ca 2+ cannot bind quickly until IP 3 has first been released, resulting in longer expected closing times. We increased the rate of IP 3 release in model 2 by a factor of 10 (that is, K 1 = 0.036 μM and K 3 = 8.0 μM), and this decreased the expected closing time. Changing the Ca 2+ concentration still has little effect on the expected closing time when the IP 3 release rate is increased.Model 3 does neither exhibit any noteworthy dependence of the cluster closing rate on the Ca 2+ concentration. Closing is obviously dominated by the Ca 2+ -independent transitions with the rate 4(l -4 + k -2 ) (see Fig.2). Fig.9A-C and Fig.10 A,Bshow that increasing the IP 3 concentration increases the expected closing time in all three models.In Fig.11A-C we show results using both a fixed and variable value for the Ca 2+ concentration, c open . When c open is fixed, it has a value of 100 μM and applies to all channels in the cluster as soon as at least one channel is open. That corresponds to very strong spatial coupling by Ca 2+ diffusion. When c open is variable it takes on different values depending on whether an individual channel is open or closed corresponding to weaker spatial coupling. For an open channel c open = 0.1 + 99.9 =

Fig. 14 A

 14 Fig. 14 A-C show how the expected closing time for a single channel varies as the number of aggregate states is reduced. When there are 330 aggregate states for models 1 and 2 and 1820 states for model 3, the full model is being used. The calculations are done for a single channel here, as it is too computationally expensive to use more channels with a large number of aggregate states. The vertical lines mark the number of aggregate states used in the results presented in this study. Note, the horizontal axes give the number of states prior to the open states being separated out. We use 12 aggregate states for model 1, 2 and 11 for model 3 here (after separating the open states). This is the maximum number for which results could be computed due

  Fig. 15 compares exact analytic results with simulations. The error does not vary greatly as the number of channels increases from 1 to 5, suggesting that if the results

A c c e p t e d m a n u s c r i p t from 20

 20 μM to 200 μM. This is due to an increased rate of Ca 2+ inhibition. In model 2 the rate for Ca 2+ inhibition while IP 3 is bound is very slow, resulting in less sensitivity of the expected closing time to the Ca 2+ concentration. Model 3 exhibits only a very weak dependence of the expected closing time on Ca 2+ . Using all three models, the expected cluster closing time increases monotonically with the IP 3 concentration.

Shuai et al. [ 12 ]

 12 estimate 40-70. However, Fraiman et al. [8] estimate 4 or 5, which agrees with an estimate by Swillens et al. [17] and Demuro et al. of about 10 channels in a cluster [20]. Very recent experiments suggest an average number of 4 [21].

A c c e p t e d m a n u s c r i p t

  Demuro, A. and I. Parker. 2008. Multi-dimensional resolution of elementary Ca 2+ signals by simultaneous multi-focal imaging, Cell Calcium 43(4):367-374 [21] Taufiq-Ur-Rahman, A. Skupin, M. Falcke, C.W. Taylor. 2009. Clustering of IP 3 receptors by IP 3 retunes their regulation by IP 3 and Ca 2+ , Nature, in press [22] Smith, F., S.M. Wiltgen, and I. Parker. 2009. Localization of puff sites adjacent to the plasma membrane: Functional and spatial characterization of Ca 2+ signaling in SH-SY5Y cells utilizing membrane-permeant caged IP 3 , Cell Calcium 45:65-76

A c c e p t e d m a n u s c r i p t

Figure 1 :

 1 Figure 1: The nine state model of an IP 3 receptor channel subunit proposed by Rüdiger et al. [10]. The Ca 2+ concentration is denoted by c and the IP 3 concentration is denoted by p. The dotted lines indicate the transitions that are made small in model 2.

Figure 2 :

 2 Figure 2: The model of an IP 3 receptor channel subunit given by Sneyed [14]. The channel states are R: unbound receptor, I 1 and I 2 : inactivated, O: open, A: activated, and S: shut. The Ca 2+ concentration is denoted by c and the IP 3 concentration is denoted by p.

Figure 3 :

 3 Figure 3: Open probability distribution F open (A, t) for various c close values. A: F open (A, t) for model 1. B: F open (A, t) for model 2. C: F open (A, t) for model 3.

Figure 4 :

 4 Figure 4: A: The expected opening times T open for the three models. B: Enlarged plot of the expected opening times for the three models to show clearly the minimum of models 1 and 2 for small values of c open .

Figure 5 :

 5 Figure 5: Opening time distribution F open (A, t) when changing p for some representative IP 3 concentrations. A: F open (A, t) for model 1. B: F open (A, t) for model 2.

Figure 6 :

 6 Figure 6: A: The expected opening times T open for the three models as function of the IP 3 concentration. B: Enlarged plot of T open for small concentrations.

Figure 7 :

 7 Figure 7: Cluster closing time distribution F close (A, t) when changing c open . A:F close (A, t) for model 1. B: F close (A, t) for model 2. C: F close (A, t) for model 3.

Figure 8 :A c c e p t e d m a n u s c r i p t

 8 Figure 8: The expected closing times T close for the three models. Model 2 is represented by the upper lines, where the dashed line with ×-symbols show the simulation results and dashed line without shows analytic results. Model 1 is shown by two

Figure 9 :

 9 Figure 9: Cluster closing time distribution F close (A, t) when changing p. A: F close (A, t) for model 1. B: F close (A, t) for model 2. C: F close (A, t) for model 3.

Figure 10 :

 10 Figure 10: A: The expected closing times T close for models 1 and 2, from both the analytic and simulation results. Model 1 is represented by the upper dashed line (analytic results) and dashed line with ×-symbols (simulation results). Model 2 is represented by the lower full lines. There are few differences between analytic (full line) and simulation results (full line with +-symbols). B: The expected closing times T close for model 3 only with analytic (full line) and simulation results (full line with +-symbols).

Figure 11 :

 11 Figure 11: F close (A, t) using both a fixed and variable Ca 2+ concentration. Here c close = 0.1 μM and p = 0.15 μM. c open = 100 μM when fixed. When c open is variable, it has a value of 100 μM for a channel which is open and a value of 10.09 μM (labeled 'Variable c open 1') or 1.099 μM (labeled 'Variable c open 2') for a channel which is closed. A: Results for model 1. B: Results for model 2. C: Results for model 3.

Figure 12 :

 12 Figure 12: Analytic and simulation results showing F open (A, t) for all three models. Here c close = 0.1 μM, and p = 0.15 μM. The full lines represent the analytic results and the -symbols the simulation results. A: Results for model 1. B: Results for model 2. C: Results for model 3.

Figure 13 :A c c e p t e d m a n u s c r i p t

 13 Figure 13: Analytic and simulation results showing F close (A, t) for the three mod-

Figure 14 :

 14 Figure 14: The effect of the number of aggregate states on the expected closing time. These results are for a single channel. A-C show the expected closing time T close as the number of aggregate states is reduced. This is shown for the three models at various parameter values. The vertical lines represent the number of aggregate states which are used in the results presented here. D-F show the number of open aggregates N open as function of the aggregate states for all three models and various parameters.

Figure 15 :

 15 Figure 15: Varying the cluster size. Graphs A, B and C show the expected closing time of the cluster T close , using models 1 (A), 2 (B) and 3 (C) respectively, as the number of channels in the cluster is varied from 1 to 5. We use p = 0.15 μM, c open = 100 μM and c close = 0.1 μM(to compute the initial condition). Both analytic (full lines) and simulation results (dashed lines with +-symbols) are shown.

Figure 2 :

 2 Figure 2:

Figure 3 :

 3 Figure 3:

Figure 5 :

 5 Figure 5:

Figure 7 :

 7 Figure 7:

Figure 9 :

 9 Figure 9:

Figure 11 :

 11 Figure 11:

Figure 12 :

 12 Figure 12:

Figure 13 :

 13 Figure 13:

Figure 14 :

 14 Figure 14:

Figure 15 :

 15 Figure 15:

A c c e p t e d m a n u s c r i p t and

  each channel state occurring at rest. Let F ch (A, t) denote the probability distribution for a channel to first open at time t given that the initial state is unknown let F open (A, t) denote the probability distribution for a cluster of channels to first open at time t (opening time distribution). Then F open (A, t), for t > 0, is given by

  in the time interval [t, t + dt] is given by the probability that one of the N

	channels is first activated in the time interval [t, t + dt] and all other channels have
	not yet been activated.
	When t = 0, if we used the method above to calculate F open (A, 0), this would
	mean that the probability for the cluster to open in the interval [0, dt] is given by the
	probability that one of the N channels opens in [0, dt], while the other N -1 do not.
	However, this is not correct. The probability that the cluster opens in [0, dt] is given
	by the probability that at least one of the N channels opens in [0, dt], while the other
	N -1 channels can either open or not open. F open (A, 0) can be calculated using

  3binding site is occupied. Specifically, if IP 3 is bound, then Ca 2+ cannot bind to the inhibitory Ca 2+ binding site, and if IP 3 is not bound, then Ca 2+ cannot bind to the activating Ca 2+ binding site. Thus, model 2 uses the same scheme as model 1, but the activating Ca 2+ binding rate when IP 3 is not bound, and the inhibitory Ca 2+

Table 1 :

 1 Parameter values for Model 1

	Parameter Value	Parameter Value
	a 0	550.0 s -1	b 0	80.0 s -1
	K 1	0.0036 μM a 1	60.0 μM -1 s -1
	K 2	16.0 μM	a 2	0.2 μM -1 s -1
	K 3	0.8 μM	a 3	5.0 μM -1 s -1
	K 4	0.072 μM	a 4	0.5 μM -1 s -1
	K 5	0.8 μM	a 5	150.0 μM -1 s -1

Table 2 :

 2 Parameter values for Model 3

	Parameter Value	Parameter Value

Several hours using a UniServer 3346 with 4Dual-Core AMD Opteron 800 series CPU's


Marco Scianna

Roeland M H Merks

Luigi Preziosi
email: luigi.preziosi@polito.it

Enzo Medico

Individual cell-based models of cell scatter of aro and mlp-29 cells in response to hepatocyte growth factor

Keywords: Cellular Potts Model, Hepatocyte Growth Factor, ARO cell line, cell-cell dissociation, invasive growth, MLP-29 cell line, scattering

The different behaviors of colonies of two cell lines, ARO (tyroid carcinoma-derived cells) and MLP-29 (Mouse Liver Progenitor cells), in response to Hepatocyte Growth Factor (HGF) are described deducing suitable Cellular Potts Models (CPM). It is shown how increased motility and decreased adhesiveness are responsible for cell-cell dissociation and tissue invasion in the ARO cells. On the other hand, it is shown that, in addition to the biological mechanisms above, it is necessary to include directional persistence in cell motility and HGF diffusion to describe the scattering and the branching processes characteristic of MLP-29 cells.

A c c e p t e d m a n u s c r i p t INDIVIDUAL CELL-BASED MODELS OF CELL SCATTER OF ARO AND MLP-29 CELLS IN RESPONSE TO HEPATOCYTE GROWTH FACTOR

Introduction and Phenomenological Description

Hepatocyte Growth Factor (HGF), originally identified as a mitogen for hepatocytes [START_REF] Mongiovi | The tyrosine kinase receptors Ron and Sea control "scattering" and morphogenesis of liver progenitor cells in vitro[END_REF][START_REF] Nakamura | Purification and characterization of a growth factor from rat platelets for mature parenchymal hepatocytes in primary cultures[END_REF], was subsequently shown to be identical to Scatter Factor (SF) [START_REF] Brinkmann | Hepatocyte growth factor/scatter factor induces a variety of tissue specific morphogenic programs in epithelial cells[END_REF][START_REF] Comoglio | HGF: a multifunctional growth factor controlling cell scattering[END_REF][START_REF] D' Alessandro | High-throughput approaches to explore the transcriptional network leading to invasive growth[END_REF], a ligand with a dramatically different activity of inducing epithelial cells dissociation [START_REF] Comoglio | Scatter factors and invasive growth[END_REF]. HGF, produced primarily by mesenchymal cells, is therefore a unique growth factor that elicits multiple cellular responses, via its receptor, a tyrosine kinase encoded by the proto-oncogene Met [START_REF] Maere | Modeling gene and genome duplication in eukaryotes[END_REF][START_REF] Matsumoto | NK4 gene therapy targeting HGF-Met and angiogenesis[END_REF].

The main biological processes regulated by HGF are mitogenesis, motility and morphogenesis [START_REF] Tamagnone | Control of invasive growth by hepatocyte growth factor (HGF) and related scatter factors[END_REF][START_REF] Trusolino | Scatter-factor and semaphorin receptors: cell signalling for invasive growth[END_REF], cell dissociation, migration through the extracellular matrix, acquisition of polarity and tubule formation [START_REF] Sakata | Heparin binding and oligomerization of hepatocyte growth factor/scatter factor isoforms[END_REF][START_REF] Stoker | Scatter factor is a fibroblastderived modulator of epithelial cell mobility[END_REF]. This combination of events, also known as "invasive growth", is fundamental during the embryonic development of most epithelial tissues. When inappropriately activated, this genetic program confers invasive ability to normal and neoplastic epithelial cells [START_REF] Mongiovi | The tyrosine kinase receptors Ron and Sea control "scattering" and morphogenesis of liver progenitor cells in vitro[END_REF][START_REF] Weidner | Scatter factor: molecular characteristics and effect on the invasiveness of epithelial cells[END_REF].

Increased levels of the HGF receptor Met have been observed in primary thyroid cells [START_REF] Lahat | Divergent effects of cytokines on human leukocyte antigen-DR antigen expression of neoplastic and nonneoplastic human thyroid cells[END_REF][START_REF] Merks | Contact-inhibited chemotactic motility: Role in de novo and sprouting blood vessel growth[END_REF], suggesting that Met may mediate some of the phenotypic alterations observed in these cells: these data allowed to precisely define the role

A c c e p t e d m a n u s c r i p t

of the HGF and its receptor in the evolution of thyroid tumors, finding which of the many cellular events, known to follow Met activation, effectively take place in thyroid cancer cells.

In this work we analyze the biological responses to HGF of a poorly differentiated thyroid carcinoma-derived cell line, ARO, that overexpresses Met [START_REF] Di Renzo | Expression of the Met/HGF receptor in normal and neoplastic human tissues[END_REF]. In these cells, HGF stimulation doesn't make the growth rate increase, probably due to the high basal level of proliferation, but it increases motility [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF].

In particular, HGF induces a mesenchymal transition with a dramatic change from a virtually non motile phenotype to a phenotype with high short range motility: cells loosen contacts, dissociate, move from their original site and start wondering in the close proximity, displaying an evident ability to invade the open spaces and modifying the interaction with the extracellular matrix (see Fig. 1 and Supp. Movie 1). Noteworthy, in response to HGF, ARO cells do not show a dramatic cytoskeletal reorganization or polarization but lose their tendency to maximize the contact area with the surrounding cells and the physiological regulation of contact-inhibition mechanism. The change in the above mentioned adhesive properties is a consequence of the change in the transmission of a number of local, contact-dependent signals (juxtracrine). In fact, we can hypothesize that the loss of intercellular adhesions does not cause by itself a dissociation phenotype, since E-cadherin -/-cells may still form strong adhesive junctions. Instead the loss of E-cadherin signaling may be responsible for this behavior [START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF][START_REF] Ramis-Conde | Modeling the influence of E-cadherin-beta-catenin pathway in cancer cell invasion: a multiscale approach[END_REF]. The above assumptions are validated also by the fact that in absence of HGF the transmission of the juxtracrine is recovered so that the dissociated, disconnected cells have the tendency to come again in contact and to form new compact colonies (not shown). Another scattering example is given by cultures of MLP-29, a strongly HGFresponsive epithelial cell line established from mouse embryonic liver. MLP-29 are small cells with ovoidal nuclei, forming tightly packed colonies.

After the addition of nanomolar concentration of HGF, the cells start to elongate, migrate and differentiate and the islands assume the characteristic stellar configuration within 1 day (as shown in Fig. 2 and Supp. Movie 2, also at at calvino.polito.it/∼preziosi/MLPex.avi). This morphogenic response is dose dependent and consists of three steps: 2. in the following 6 h there is a cytoskeletal reorganization, followed by an asymmetric polarization and an active migration of cells, which lead to the formation of some sprouts;

3. within 24 hours, at the end of the branching process, the colonies assume the stellar shape, with sprouts formed by one or few elongated cells that keep a thin adhesion contact with the rest of the island.

None of these events is observed in control medium or after stimulation with other growth factors: in those cases, MLP-29 cells tend to form spherical and compact aggregates of undifferentiated, almost round cells, with the interstitial spaces almost filled with microvilli protruding from the cell surface along the entire perimeter of each cell [START_REF] Mongiovi | The tyrosine kinase receptors Ron and Sea control "scattering" and morphogenesis of liver progenitor cells in vitro[END_REF].

Entering more in details in the description of HGF-driven scattering, the cells in the middle of an island don't have a dramatic change of shape, keeping their original elliptical shape, and passively move only dragged by the nearest sprout.

Instead, the cells at the margin acquire a polar configuration driven by an internal remodeling of the actin cytoskeleton: they are highly asymmetric and have distinct apical, lateral and basal surfaces which correspond to a leading and trailing surface during their migration. In particular, they migrate faster along their long axes than along their short axes, because they need fewer surface protrusions and retractions to move normal to their short interface than to their long interface. Such anisotropic, cell-shape-dependent motility results in strong

A c c e p t e d m a n u s c r i p t

directional persistence because the cells must reorient before moving into a new direction. This inertial dynamic introduces two timescales into scattering: fast migration along the cells' long axes produces rapid branching and polarization, the slower sideways migration coarsens the pattern as the cells align and close small lacunae. The polarized cells, in the part of their surface in contact with the ECM, extend filopodia towards the open space and keep a thin adhesion point with the main corpus of the colony. The branches, generally formed by couples of parallel elongated cells, are radially disposed with angles ranging between 15 and 100 degrees. If the HGF is no longer added in the culture a sort of involution of the sprouts is seen and the colonies tend to become round and compact again (not shown).

We here model the just described behavior of both ARO and MLP-29 cell lines at a mesoscopic level using the Cellular Potts Model (a lattice-based, Monte Carlo approach developed for instance in [START_REF] Graner | Simulation of biological cell sorting using a twodimensional extended Potts model[END_REF][START_REF] Merks | Dynamic mechanisms of blood vessel growth[END_REF][START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF][START_REF] Merks | Contact-inhibited chemotactic motility: Role in de novo and sprouting blood vessel growth[END_REF][START_REF] Merks | The Glazier-Graner-Hogeweg model: extensions, future directions, and opportunities for further study[END_REF]), which allows to describe cell-cell adhesion and cell migration. An energy minimization philosophy, a set of constraints and auxiliary conditions determine how the cells move. In particular intercellular junctions determine an adhesive (or binding) energy between cells and an energy constraint regulates cell volume.

We show also how, with respect to ARO colonies, in order to describe the behavior of MLP-29 colonies it is necessary to include in the model a specific term mimicking cell persistence. In both cases, we assume that the HGF is not limiting and distributes uniformly over the ECM.

The results of the simulations are in good agreement with the experiments obtaining all the described features and morphologies characterizing the evolution of both cell lines .

Cellular Potts Model

We use the CPM and make all of our biological assumptions at the cell level, while discussing the model results at the level of the tissue. Our basic model is referred to the case of the ARO cells and assumes that (a) the cells keep a round shape, (b) the addition of the growth factor influences the cellcell adhesion energy and the system motility, leading to the dissociation and invasion processes and (c) the E-cadherin binding at cell-cell interfaces represses the formation of surface projections (filopodia), while they form normally at the unbound parts of the cell surface in contact with the extracellular matrix (contact-inhibition).

The CPM represents biological cells as patches of lattice spins σ(x) of identical state on a second-nearest-neighbor square lattice, where each state labels a single biological cells, see Fig. 3A. Connections between neighboring lattice sites of unlike state σ(x) = σ(x´) represent membrane bonds between cells, where the non-dimensional parameters J σ(x),σ(x´) give the bond energy, assuming that the types and the number of E-cadherins and other adhesive cell-surface proteins determine the Js. Since in our model there is only a type of cells placed in a virtual extracellular matrix the Js reduce to J c,M and J c,c , where c represents a cell and M the ECM, Fig. 3B. We further define the relative adhesion through the parameter γ c,M = J c,M -J c,c /2, which enables us to determine whether the cells cohere (γ c,M < 0) or dissociate (γ c,M > 0). An energy penalty increasing with the cells' deviation from a designed target volume, A σ , constraints the volume of the cells a σ , i.e. A σ and a σ are expressed in lattice sites. We express the adhesion and target area deviations in terms of the effective energy H:

A c c e p t e d m a n u s c r i p t

H = x,x´J σ(x),σ(x´) (1 -δ σ(x),σ(x´) ) + λ σ (a σ -A σ ) 2 , ( 1
)
where x´represents the eight second-order neighbors of x, the non-dimensional Lagrangian multiplier λ the cells' resistance to compression and the Kroenecker delta is δ x,y = {1, x = y; 0, x = y}. The cells reside in an ECM which we model as generalized CPM cells without a volume constraint and with σ = 0. We define a special, high cell-border energy J c,B = 100 to prevent cells from adhering to the boundaries of the simulated area (which we define as frozen pixels of type B).

To mimic cytoskeletally driven, active surface fluctuations, we randomly choose a lattice site, x, and attempt to copy its spin σ(x) into a randomly chosen neighboring lattice site x´. For better isotropy we use the twenty, first to four nearest neighbors. This attempt represents a lamellipodial or filopodial protrusion of the first cell displacing the surface of the second cell. During a Monte Carlo Step (MCS), which can be evaluated to be 30 seconds long from a direct comparison of simulations and experiments, we attempt a number of copies equal to the number of sites in the lattice. We minimize H using standard Metropolis dynamics: we calculate how much the free energy H would change if we performed the copy and accept the attempt with Boltzmann probability:

P (ΔH) = ⎧ ⎨ ⎩ e -ΔH-H 0 T ΔH ≥ -H 0 1 Δ H < -H 0 (2)
where H 0 is an energy threshold which models viscous dissipation and energy loss during bond breakage and formation and the non-dimensional coefficient T is related to cell motility. Each CPM cell i moves according to the effectiveenergy gradient around it, v i ∝ ∇ i H, where ∇ i H is the component of the gradient relative to cell i. Thus, given an effective energy, we can calculate the resulting cell motion and the force required to create such a motion. In the

A c c e p t e d m a n u s c r i p t

highly viscous environment of the ECM, the drag force F i exerted on a cell i is proportional to its velocity v i : this relation is often called the overdamped force-velocity response. Thus the force exerted on a cell is proportional to the effective energy gradient

F i ∝ ∇ i H.
In the case of ARO aggregates, HGF stimulation is modeled only through its effects on cell-cell adhesion energy and on cellular motility. In fact we can preliminary hypothesize a directed proportionality between the virtual HGF concentration and both the model parameters T and J c,c , so that increasing doses of the growth factor correspond to increasing values of cellular motility T and of inter-cellular adhesion coefficient J c,c , mimicking weaker cell-cell junctions. J c,M remains fixed in each simulation, since in this model the matrix is only considered as a rigid thin layer formed by generalized cells, that are not represented in the CPM lattice and not influenced by HGF.

As we shall see, the basic model ( 1) is not able to predict the behavior of the MLP-29 cells. In fact, in a culture of MLP-29 cells, the growth factor induces also cell elongation and polarization. These ingredients, that are not contained in [START_REF] Brinkmann | Hepatocyte growth factor/scatter factor induces a variety of tissue specific morphogenic programs in epithelial cells[END_REF], are essential for branching morphogenesis and the formation of the stellar shape aggregates experimentally found.

In fact, still trying to keep contact with the colony, MLP-29 cells take time to reorganize their internal machinery to change their direction of motion and their orientation and polarization induce a persistence motion, as done in [START_REF] Glazier | The Glazier -Graner -Hogeweg Model: Extension, Future Direction, and Opportunitites for Further Study. Single-Cell-Based Models in Biology and Medicine[END_REF]. We can model this cytoskeletal persistence time by adding to our effective energy H an inertial-like term of the form:

H in = λ in x,x´| |v(σ, t) -v(σ, t -Δt)|| 2 (3)
where v(σ, t)

v(σ, t) = cm(σ, t) -cm(σ, t -Δt) Δt , ( 4
)
is the instantaneous center of mass (cm(σ) = a -1 (σ) x:σ(x)=σ x) velocity of cell σ and the non-dimensional Lagrangian multiplier λ in controls the persistence time. Typically, Δt is one or more MCS. Obviously, if λ in = 0 cells undergo uncorrelated Brownian motion, while if λ in is large the motion is almost ballistic.

In this case, it is also needed to add a PDE to model the diffusion and absorption of HGF in the culture from a macroscopical point of view [START_REF] Merks | Contact-inhibited chemotactic motility: Role in de novo and sprouting blood vessel growth[END_REF][START_REF] Merks | The Glazier-Graner-Hogeweg model: extensions, future directions, and opportunities for further study[END_REF]: the equation regulating the evolution in time of the growth factor is

∂h ∂t = α -εδ σ(x),0 h + D∇ 2 h ( 5
)
where

δ σ(x),0 = ⎧ ⎨ ⎩ 1 inside the cell 0 in the ECM, (6)
α is the rate at which the HGF is uniformly added, ε is the uptake rate of the growth factor and D the diffusion constant, which will be always set to 10 -13 m 2 s -1 , as done in [START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF]. In every experiment, the secretion and the decay

A c c e p t e d m a n u s c r i p t

rate are set equal, see also [START_REF] Merks | Contact-inhibited chemotactic motility: Role in de novo and sprouting blood vessel growth[END_REF]: for our parameter settings the HGF diffuses more rapidly than the cells move, enabling us to ignore advection which may occur as the cells displace the ECM. We implement chemotaxis due to preferential extension of pseudopods up HGF gradients by including a chemical energy change at each copy [START_REF] Merks | Dynamic mechanisms of blood vessel growth[END_REF][START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF][START_REF] Merks | Contact-inhibited chemotactic motility: Role in de novo and sprouting blood vessel growth[END_REF]:

ΔH chemotaxis = -μ (h(x´) -h(x)) , ( 7
)
where h is the concentration of HGF, which, as said, we assume is present everywhere in the layer of HGF under the cells, x´is the target site and x the source site. In this case chemotaxis depends on the HGF gradient only, independently of its concentration.

The chemotaxis coefficient μ is set equal to 0 at cell-cell interfaces and to 500 at cell-ECM interfaces, ensuring that chemotactic extensions occur only at cell-ECM surfaces, mimicking contact-inhibition chemotaxis.

The PDE is solved numerically using a finite difference scheme on a lattice that matches the CPM lattice, using 15 diffusion steps per MCS, i.e. the time step for the diffusion solver is 2 sec and the space step is 2 μm.

To sum up, the scheme used for MLP-29 cell line consists of an hybrid CPM/PDE model: the CPM models the cells, the PDE models the diffusion and the absorption of the growth factor.

Also in the MLP case, each Monte Carlo Step (MCS) corresponds to 30 s: for this choice of timescale, the mean of the cell velocities agrees well with in vivo observations.

Simulation of the Scattering Process of ARO Cells

Figure 4A-D shows a typical simulation of the time-sequence evolution of a culture of ARO cells in presence of HGF: we initiate the simulation with a mass of 16 virtual cells over an area of 300 × 300 lattice sites, which we positioned in a larger lattice of 500 × 500 to minimize boundary effects. Each CPM lattice site corresponds to a square of size 2 μm × 2 μm. We assume that cells do not divide or grow during patterning and set the target area of the simulated ARO cells to 50 lattice site, corresponding to 200 μm 2 .

Based on cell culture measurements [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF][START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF], we set λ = 25, J c,c = 60, J c,M = 20 and T =50. After an initial roughening, the cells start to detach and to radially spread in close proximity, keeping an almost elliptical shape. In Fig. 4E-H we also simulate the loss of HGF in the culture medium by decreasing J c,c to 20, maintaining the other parameters fixed: as in the experiments, the cells aggregate again in an almost round compact island, confirming the hypothesis of the recovered E-cadherin signalling.

The evolution of the simulated colony coincides well with the experimental evidences, having also the same temporal dynamics, as shown in Fig. 5, where we plot the Scattering Index

S I (t) = A(t) A(0) -1, (8)
with A(t) defined as the area of the minimal convex polygon with vertices at the more external cells (at the tips of the branches in MLP29 case) at time t and A(0) at t = 0. The experimental scatter begins slightly slower than the simulated, however, after the initial phase (t > 1500 MCS), S I increases at comparable rates in both cases.

A c c e p t e d m a n u s c r i p t

To further study the robustness of the model and to identify the critical features of the output components, we analyze how the results of the simulations depend on the parameters by measuring the Scattering Index after 3000 MCS (i.e., 1 day). Fig. 6 shows that, at suitable T , our simulated aggregates do not scatter for low values of inter-cellular adhesion parameter J c,c , corresponding to strong inter-cellular adhesion, whereas they dissociates when cells adhere less strongly: this observation corresponds well with our experiments where ARO aggregates scatter at high concentration of HGF, see [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF]. In Fig. 6 we also look for the motilities required for the process: for T < 30 there is no detectable scattering within the time limit (3000 MCS), independently from the adhesion energies J(c, c) (the cells dissociate, but don't spread away), for T > 60

A c c e p t e d m a n u s c r i p t

they break up in small pieces or fragments, a well-characterized non-biological artifact of the CPM [START_REF] Merks | The Glazier-Graner-Hogeweg model: extensions, future directions, and opportunities for further study[END_REF], while for the in-range values, the process becomes more significant and quicker for every T-increment. These results validate our hypothesis that HGF concentration can be modeled increasing simultaneously T and J(c, c), since both an increasing in cell motility and the brokage in E-cadherin bonds are needed for ARO cells scattering, as seen in experimental evidences.

Figure 7: HGF-mediated increase of motility: comparison of single-cell motion tracks obtained by 24-h time-lapse video-recordings (A-B, see [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF]) and 3000 MCS simulations (C-D) of ARO cells cultured in the absence (N.S.) and in presence of HGF, as indicated. In the experiments the trajectories are obtained by tracking continuously the motion of the nucleus of a cell, while in the simulations, they are formed by the position of the center of mass of the single cells at each time step. Experimentally, this would correspond to connect the position of the nucleus every 30 seconds (i.e., 1 MCS is set equal to 30 sec). Ten representative tracks were chosen for each case and optimally arrayed for picture presentation. The oval positioned at the center of each group of tracks represent the size of a single cell. Fig. 7 shows a definitive confirmation of the HGF-induced increase of cell motility, comparing the single-cell trajectories during a 24 h (= 3000 MCS) time-lapse, both in the control and in the stimulated case resulting from the experiment and the numerical simulation: the value of the mean effective distance covered by an ARO cell is about 5 μm in the absence and 12 μm in the presence of the growth factor (means over ten representative cells, error of ≈ 1 μm in both cases).

These evidences are used also to present the basic features of a simulation reproducing a "wound healing" assay, which tests the ability of a cell line to fill gaps created in cell monolayers and is generally considered a simple and reliable test for quantitative evaluation of cell motility. Experimentally, a confluent monolayer of ARO cells, grown in fetal calf serum, is wounded with a pipette tip and incubated for 24 h in the presence HGF, see [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF]. In the simulation, Fig. 8, the culture is formed by two masses of 120 virtual cells distributed over an area of about 700 μm × 700 μm (350 × 350 lattice sites), while the other parameters are the same as in the simulations of Fig 4 . The gap is of about 120 μm (60 lattice sites): the maximum value, deduced from the above results, for which the two colonies encounter within the time limit (24 h). The time needed to invade the space between the colonies is evaluated in Fig. 9

A c c e p t e d m a n u s c r i p t

Simulation of the Scattering Process of MLP-29 Cells

Figure 10 shows a typical simulation of a time-sequence evolution of a colony of MLP-29 cells in the presence of HGF: we initiate the simulation with a mass of 16 virtual cells, each of them with an area of ≈ 160 μm 2 (40 lattice sites), over an area of ≈ 700 μm ×700 μm (333 × 333 sites), which is positioned inside a larger lattice of 1000 μm × 1000 μm (500 × 500 sites) to minimize the boundary effects. Assuming that the cells do not divide during patterning and in accordance with culture measurements [START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF], for this sample simulation the HGF's secretion and decay rate are set equal to 10 -4 s -1 . In this and in the other simulations we set the effective cell-cell adhesion equal to the cell-ECM adhesion (i.e. J c,c = 2J c,M , the factor of 2 arise because we model ECM as a single large generalized cell), since experimental analysis suggests a reorganization of intercellular bonds rather than a completely brokage of them. The other parameters are set, in this case, equal to λ = λ in = 20 and T = 50.

After an initial stage in which the surface of the colony roughens, with some cells radially protruding from the surface, some branches develop still keeping contact with the neighboring cell and the center of the island: this gives the colony the stellar-shape configuration experimentally found.

As soon as the HGF is no longer added (after 24 hours ≈ 3000 MCS), the sprouts start a sort of involution and the island returns compact, a biological behavior that is possible to see also experimentally in [START_REF] Mongiovi | The tyrosine kinase receptors Ron and Sea control "scattering" and morphogenesis of liver progenitor cells in vitro[END_REF], see Fig. 10. We model the decay of HGF by imposing α = 0 in Eq. ( 5) (before α = 10 -4 s -1 ), while the other parameters are kept unchanged. The evolution of the simulated colony coincides well with the experimental evidences not only qualitatively as shown in Fig. 10, but also quantitatively as it is possible to see in Fig. 11, that compares the time evolution of the Scattering Index of the island as defined in ( 8) for the experiment and the simulation: it is worth noticing that also the simulated branching develops through the three temporal phases experimentally pointed out in the introduction. The Scattering Index does not significantly increase after ≈ 2500 MCS since it is possible to hypothesize an equilibrium between the different forces: the tendency of the cells to elongate and to respond to HGF can no more overcame the adhesion with the rest of the colony and the area constraint. A breakage in this equilibrium would make the scattered island fragment. An analogous behavior is experimentally seen, with the scattered colonies that reach and keep a limit conformation before the detachment of some branches. The simulations, confirming the experimental evidences, show a sort of differentiation of the cells in the colonies: those in the sprouts, after the polarization, have a persistence motus in the direction of the branch they belong to, the others (in the valley or in the center) keep a round shape and have a more isotropic To have a deeper idea of the quantity of HGF needed for sprouting and to test whether its dose affects the results, the simulation model is run with increasing values of α (= ε, i.e.), while the other parameters are unchanged. The initial island doesn't sprout for α < 10 -4 s -1 , requiring a minimal level of HGF to start the process, while a sort of phase transition is observed (see Fig. 13) for α > 10 -3 s -1 : between these two values branching process regularly develops.

The tendency of the cells to polarize, elongate, and have a persistent dynamics affects the geometry of the branches. This is studied running a set of simulations keeping the same initial configuration as that in Fig. 10 and gradually increasing λ in (the other parameters are unchanged, with α = ε = 10 -4 s -1 ). For 1 < λ in < 45 the sprouts became longer and thinner with an almost linear trend respect to the λ in -increments. In fact for very low values of the parameter they are almost formed by four or five parallel cells which are not completely polarized, Fig. 14A-B. For λ in > 45 the branches tend to become only one cell wide and constituted of two or three elongated cells in length. For still higher values, there is a sort of limit threshold: obviously the branches can't be formed by less than a single cell and those which are longer than ≈ 52 μm tend to detach from the main corpus of the colony and, as isolated cells, to migrate in the direction of their longer axes, a behavior also observed in the experiments.

It is also worth noticing that the geometry of the branches is not dramatically influenced by the dose of HGF, since the curves related to different quantities of the growth factor have almost the same slope and don't differ much in the particular values, but may be also related to the number of cells forming the studied island (as shown in Fig. 14C). For this reason some simulations are performed starting with colonies with different densities.

As the number of cells increases, the lattice is taken wider, to avoid overlapping problems. For low densities (n < 60) the branches are quite long and thin, so that it can be hypothesized that the colony assumes a stellar-shape configuration, with a little main corpus. For high densities (n > 60) the sprouts tend to become shorter and thicker, while at still higher densities (n > 130) there are limit thresholds for both length and thickness. The decrement of the length and the relative increment of the thickness is almost specular and has the consequence that the area of a virtual branch remains almost constant, independently from the cell density, which can only lead to a bigger main corpus of the island. The cell density affects also the number of the branches developing in a timelapse of 24 h: for very low values (n < 8) there is no development of sprouts, for high densities (n > 108) there is a constant value for their quantity, as shown in Fig. 14D. In the medium range, the number of the branches increases with the number of cells forming the island. The higher quantity of sprouts appearing in bigger colonies may have two explanations: first, bigger colonies have more cells available for branching, then they have wider external surfaces and so more space for the growth of the sprouts.

In our model we have assumed that the MLP29s' decision to extend or retract a pseudopod depends on the concentration difference between the retracted and extended position of the pseudopods, regardless of the absolute concentration of the HGF. However, at higher concentration we can hypothesize that most HGF receptor will saturate, thus becoming insensitive to HGF levels, a situation that occurs experimentally and that is possible to find also in the sprouting process of blood vessel, [START_REF] Merks | Dynamic mechanisms of blood vessel growth[END_REF][START_REF] Merks | Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling[END_REF]

A c c e p t e d m a n u s c r i p t Discussion

We have studied the scattering processes induced by HGF on both ARO and MLP-29 cell lines and shown that the related simulations are in good agreement with experimental results. We have demonstrated that a change in cell-cell adhesion properties, in conjunction with an increment in the system motility, suffices to explain aspects of the biological responses triggered by HGF on a culture of ARO cells. In the absence of HGF the cells are characterized by strong cell-cell junctions, as they form closely adjoined layer, and by a virtual non-motile phenotype: in particular, under normal conditions, ARO cells don't usually detach and move away from the nearest neighbors. On the contrary, HGF induces both the dissociation and invasion programs. This might be related to an interplay with the activity of the cell-adhesion molecule E-cadherin and with its downstream signaling pathways, e.g. via BCL9-2, [START_REF] Brembeck | Essential role of BCL9-2 in the switch between beta-catenin's adhesive and transcriptional functions[END_REF].

In particular, the results of the simulations show that round ARO cells isotropically attract surrounding cells and form rounded, disconnected islands but they haven't acquired autonomous patterning ability (they don't self-organize into polygonal patterns).

The distinctive effects of HGF on ARO cells, i.e., induction or enhancement of motility and invasive growth, are clearly detected at a certain range of doses, since they are not elicited for lower quantities. Therefore it can be concluded that the ARO cells acquire not only an HGF-dependent phenotype but an HGF dose-dependent one.

In the case of a culture of MLP-29 cells, we have shown that the basic model, with a change in cells' adhesive properties and an increased motility, is not sufficient to explain their behavior. It is necessary to add a term related to cell polarization and persistence, in conjunction with a chemotactic response, which were shown to be key ingredients also in the formation of vascular networks [START_REF] Gamba | Percolation, morphogenesis, and Burgers dynamics in blood vessel formation[END_REF][START_REF] Serini | Modelling the early stages of vascular network assembly[END_REF][START_REF] Tosin | Mechanics and chemotaxis in the morphogenesis of vascular networks[END_REF].

Mimicking contact-inhibition of chemotaxis by suppressing chemotaxis at cell-cell interfaces, our results suggest that the presence of the growth factor produces a continuous, inward, normal force at the cluster boundary, creating a buckling-like instability and that the chemotactic forces compress small initial bumps laterally, producing sprouts. In fact MLP-29 cells are more likely to extend pseudopods than cells in the valley between the bumps.

In other words, cells in growing tips are also more influenced by HGF than those in valleys between tips and in the center of the island: filopodia at growing tips are more frequent because they have a lower energy cost, making the rate of pseudopod extension critical to pattern evolution. For the instability to persist it is required that the cells in the valley must retract, while the others protrude, making the branch they belong to develop.

At equilibrium the tracks of the zones of HGF-action have a quasi-Gaussian profile. It levels off towards the clusters' boundary, while its inflection point is at the cell-cluster boundary.

The simulations suggest that the addition of the growth factor influences the general motility of the culture. The dissociation and the migration of isolated cells permits to study the dynamic of a generic MLP-29. The anisotropic migration requires no extra assumptions in the CPM and produces long persistence time in the presence of HGF cues, because, as seen, cells change direction slowly. The long persistence time introduces two time scales into the branching mor- The model has also helped to find the role of the cell shape in the geometry of the branches: the more a cell is elongated and tends to have a persistence dynamics in motion, the more the branch it belongs to is long and thin. Continuum and cell oriented approaches show a correlation between the density of the culture and the equilibrium geometry of the branches: thicker and shorter sprouts form in cell cultures with higher numbers of cells.

In our ongoing work we are refining and validating our computational model by obtaining experimentally derived values for the model parameters and by comparing our simulation results quantitatively to time-lapse videomicroscopy experiments. We also plan to model the interaction of the individual cells with the possibly viscoelastic matrix, which experiments and theory suggest to be crucial in aspects of tissue invasion. An alternative model that can be extended to cover this aspect is presented in [START_REF] Chauviere | A continuum model for mesenchymal motion in a fibrous network[END_REF][START_REF] Chauviere | Modelling the motion of a cell population in the extracellular matrix[END_REF].

In particular it is possible to discuss how its stiffness and elasticity can explain the relation between the HGF dose and the culture: it is possible to assume that ECM may trap the growth factors to generate a delayed scattering with unstable branches.

From the experimental viewpoint, further experiments are needed to compare morphologies in experiments using different doses of HGF.

Another interesting development is to analyze the dynamics of scattering in other cell lines under different physiological and pathological conditions and compare them to the computational model's dynamics for similar simulated conditions to find a single Hamiltonian which regulates the behavior of different cell lines in response to HGF.

Figure 1 :

 1 Figure 1: Phase contrast time-lapse microscopy of in vitro dissociation and migration of ARO cells, grown in 10 per cent fetal calf serum and incubated in the presence of 100 U/ml of hepatocyte growth factor for 2 h, 12 h, 18 h and 24 h. Scale bar, 25 μm long. The full time-lapse movie is provided in supplemental material and available in calvino.polito.it/∼preziosi/AROex.avi.

A c c e p t e d m a n u s c r i p tFigure 2 :

 2 Figure 2: Typical time sequence of in vitro branching morphogenesis of a colony of MLP-29 cells, incubated for 24 h in DMEM -1 % FCS, containing recombinant HGF, diluted to the final concentration of 100 ng/ml, at 1 h, 6 h, 10 h, 15 h, 20 h and 24 h after incubation. Scale bar is 25 μm long. The full time-lapse movie is provided in supplemental material and available in calvino.polito.it/∼preziosi/MLPex.avi

Figure 3 :

 3 Figure 3: (A) Typical square grid of the CPM. The single cells are labeled by the same spin. (B) Membrane bonds between different cells (color lines). In each picture the ECM is represented by the lattice sites of index M.

Figure 4 :

 4 Figure 4: Cell-cell dissociation process in a simulation initiated with a clump of 16 ARO cells: (A) 6 hours, (B) 12 hours, (C) 18 hours, (D) 24 hours, (E) 30 hours, (F) 36 hours, (G) 42 hours and (H) 48 hours. After 3000 MCS (≈ 24 hours) the HGF is no longer added in the culture. Full simulation is available at calvino.polito.it/∼preziosi/AROsim.avi.

Figure 5 :

 5 Figure 5: Comparison of the evolution in time of the Scattering Index of ARO cells for experimental (black dotted line) and simulation (red dotted line) evidences. Error bars show standard deviations over ten representative colonies for the experiments and over ten simulations.

A c c e p t e d m a n u s c r i p tFigure 6 :

 6 Figure 6: Dissociation process in response of HGF in a culture of ARO cells: Scattering Index at t=3000 MCS, defined as (8), as a function of (A) both cell-cell adhesion energy Jc,c and cellular motility T ; (B) only of Jc,c at T = 20 (black solid line), T = 30 (black dashed line), T = 40 (black dotted line), T = 50 (black dash-dotted line) and T = 60 (black dot-dot-dashed line); (C) of T at Jc,c = 10 (black solid line), Jc,c = 30 (black dashed line), Jc,c = 50 (black dotted line) and Jc,c = 70 (black dash-dotted line). (D) Evolution in time of the Scattering Index for T = 20 (black solid line), T = 50 (black dashed line), T = 60 (black dotted line), T = 80 (black dash-dotted line). All standard deviations over ten simulations are about 6 per cent of the relative mean values.

  Fig.7shows a definitive confirmation of the HGF-induced increase of cell motility, comparing the single-cell trajectories during a 24 h (= 3000 MCS) time-lapse, both in the control and in the stimulated case resulting from the experiment and the numerical simulation: the value of the mean effective distance covered by an ARO cell is about 5 μm in the absence and 12 μm in the presence of the growth factor (means over ten representative cells, error of ≈ 1 μm in both cases).These evidences are used also to present the basic features of a simulation reproducing a "wound healing" assay, which tests the ability of a cell line to fill gaps created in cell monolayers and is generally considered a simple and reliable test for quantitative evaluation of cell motility. Experimentally, a confluent monolayer of ARO cells, grown in fetal calf serum, is wounded with a pipette tip and incubated for 24 h in the presence HGF, see[START_REF] De Luca | Met overexpression confers HGF-dependent invasive phenotype to human thyroid carcinoma cells in vitro[END_REF]. In the simulation, Fig.8, the culture is formed by two masses of 120 virtual cells distributed over an area of about 700 μm × 700 μm (350 × 350 lattice sites), while the other parameters are the same as in the simulations of Fig 4.The gap is of about 120 μm (60 lattice sites): the maximum value, deduced from the above results, for which the two colonies encounter within the time limit (24 h). The time needed to invade the space between the colonies is evaluated in Fig.9as a function of different cell-cell adhesion parameter values, in the cases of T = 40 and T = 50. For strong inter-cellular bonds, corresponding to J c,c < 30, ARO cells display

Figure 8 :

 8 Figure 8: "Wound-healing" process: evolution in time of a simulation initiated with two masses of 120 ARO cells: (A) 1 hour,(B) 8 hours, (C) 12 hours, (D) 16 hours, (E) 20 hours and (F) 24 hours. Full simulation is available at calvino.polito.it/∼preziosi/AROsim2.avi.

Figure 9 :

 9 Figure 9: Time for wound healing as function of the cell-cell adhesion energy, for T=40 (red dotted line) and T=50 (black dotted line). Error bars show standard deviations over ten simulations.

Figure 10 :

 10 Figure 10: Branching morphogenesis in a simulation initiated with a clump of 16 MLP-29 cells. After 3000 MCS ≈ 24 hours the HGF is no longer added in the culture and regression of cell scattering is observed. Isolines (green) indicate ten HGF levels relative to the maximum concentration in the simulation. Full simulation is available at calvino.polito.it/∼preziosi/MLPsim.avi.

Figure 11 :

 11 Figure 11: Comparison of the time-evolution of the Scattering Index of a round colony of 16 cells for experimental (red dotted line) and simulation (black dotted line). The parameters of the simulations are those of Fig. 10. Error bar show standard deviations over ten representative islands for the experiments and over ten simulations.

Figure 12 :

 12 Figure 12: Hepatocyte growth factor-mediated increase of motility as shown by manual time lapse analysis. Single cell motion tracks obtained from 24 hours-lapse simulations of MLP-29 cells culture in the absence (N.S.) and in the presence of HGF. The red oval positioned at the center of the group of tracks represent size of a single cell.

Figure 13 :

 13 Figure 13: (A) Scattering Index of 16-MLP29 cell clusters after 3000 MCS (≈ 24 h) as a function of HGF secretion rate α. Error bars show standard deviations over ten simulations. Insets are representative configurations after 3000 MCS. (B) Comparison of time evolution of S I for α = 10 -5 s -1 (black solid line), α = 10 -4 s -1 (black dashed line), 2 • 10 -4 s -1 (black dotted line), 5 • 10 -4 s -1 (black dot-dashed line) and 10 -3 s -1 (black dot-dot-dashed line). All standard deviations over ten simulations are about 5 per cent of the relative mean values.

A c c e p t e d m a n u s c r i p tFigure 14 :

 14 Figure 14: Mean length (A) and mean thickness (B) of the branches formed in a 16 MLP-29 cells round island after 3000 MCS (≈ 24 h) as a function of the persistence penalty coefficient λ in in the cases of α = ε = 10 -4 s -1 (black dashed lines), α = ε = 2 • 10 -4 s -1 (black dot-dashed lines) and α = ε = 5 • 10 -4 s -1 (black dotted lines). All standard deviations over ten simulations are about 4 per cent of the relative mean values. (C) Mean length (red dotted line) and mean thickness (black dotted line) as functions of the initial number of cells forming the colony for α = ε = 2 • 10 -4 s -1 . (D) Number of the branches formed in a time-lapse of 3000 MCS (≈ 24 h) as a function of the initial cell density of the island. Error bars indicate standard deviations over ten simulations. Insets are representative configurations after 3000 MCS.

A c c e p t e d m a n u s c r i p t

  phogenesis: the faster along the longer axes produces the growth of the spikes, the slower, radially directed, coarsens the main corpus of the island.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t Figure 1 A c c e p t e d m a n u s c r i p t Figure 2 A c c e p t e d m a n u s c r i p t Figure 3 A c c e p t e d m a n u s c r i p t Figure 4 A c c e p t e d m a n u s c r i p t Figure 5 A c c e p t e d m a n u s c r i p t Figure 6 A c c e p t e d m a n u s c r i p t Figure 7 A c c e p t e d m a n u s c r i p t Figure 8 A c c e p t e d m a n u s c r i p t Figure 9 A c c e p t e d m a n u s c r i p t Figure 10 A c c e p t e d m a n u s c r i p t Figure 11 A c c e p t e d m a n u s c r i p t Figure 12 A c c e p t e d m a n u s c r i p t Figure 13 A c c e p t e d m a n u s c r i p t Figure 14 A c c e p t e d m a n u s c r i p t Figure 15

 123456789101112131415 Figure 1

A c c e p t e d m a n u s c r i p t