Simulation of melanoblast displacements reveals new features of developmental migration

ABSTRACT To distribute and establish the melanocyte lineage throughout the skin and other developing organs, melanoblasts undergo several rounds of proliferation, accompanied by migration through complex environments and differentiation. Melanoblast migration requires interaction with extracellular matrix of the epidermal basement membrane and with surrounding keratinocytes in the developing skin. Migration has been characterized by measuring speed, trajectory and directionality of movement, but there are many unanswered questions about what motivates and defines melanoblast migration. Here, we have established a general mathematical model to simulate the movement of melanoblasts in the epidermis based on biological data, assumptions and hypotheses. Comparisons between experimental data and computer simulations reinforce some biological assumptions, and suggest new ideas for how melanoblasts and keratinocytes might influence each other during development. For example, it appears that melanoblasts instruct each other to allow a homogeneous distribution in the tissue and that keratinocytes may attract melanoblasts until one is stably attached to them. Our model reveals new features of how melanoblasts move and, in particular, suggest that melanoblasts leave a repulsive trail behind them as they move through the skin.


Potential functions
We assign to each Kc a "potential" ( ) at time , which can be interpreted as a factor density but which can also have an abstract mathematical meaning of the capacity of attraction of Mb by Kc. For a free-Kc, we assume that the time evolution of Ψ is modeled by the differential equation These Kc potentials enable us to define a density of factors by a diffusion process of the factor. At each time step the Kc potential diffusion is obtained by repeating (up to three times) the following process: -At each polygon vertex, a potential is calculated by an average of the potential of the Kc to which the vertices belong. -The new Kc potential is the average of the potential of all its vertices.
Thus the factor is propagated over three Kc away after which a value for the potential is computed as a weighted average of the initial value and of the new one. If the weight ! is 0, there is no diffusion and if it is 1, the diffusion is maximal. This density of factors enables to calculate a potential of attraction ! t at each node as an average of potentials of Kc lying around the node . The potential of repulsion between Mb, denoted ! , computed at a given node , is equal to 0 or 1 according to the fact that all nodes adjacent to are occupied by another Mb or not. To model the additional assumption that a Mb can repulse another Mb for a short time after it has left a node, we consider more generally that at a given time, a node occupied by a Mb k time steps before this tim has a potential ! where 0 < < 1. Therefore, the potential ! is computed at a node as the minimum of the potential of all nodes adjacent to (with the possibility of considering several layers of adjacent nodes). Some of the experimental sets of trajectories we want to simulate exhibit a clear anisotropy of the trajectory directions without a specific orientation along this direction, possibly due to some experimental and non-biological cause. To take this anisotropy into account, a term ! ( ) is added to the potential in order to favour one direction. This term is proportional to the absolute value of cos , where is the angle between the displacement direction and the favoured direction.

Melanoblast velocity
The mean speed of a Mb is the trajectory length divided by the time duration, and in a simulation only the number of time steps is known but not the duration of a time step. Thus, we determine the time step duration in order that the average of all Mb mean speeds is equal to a given experimental value. Let us remark that if the number of time steps is large enough, all of the trajectory lengths will be almost equal to , where is the mean value of the Kc polygon side lengths. In this case, the distribution of lengths is very narrow and only depends on geometric properties (e.g. the Kc polygonal shapes). However, the trajectory lengths must also depend on biological features of the Mb. The classical way to correct this model drawback is to allow for the Mb sometimes to be stationary. That is why in order to obtain a realistic distribution of Mb mean speeds the following process is added: a Mb numbered is assigned a probability of moving ! with a set of probabilities ! uniformly distributed between a value !"# and 1. At each time step and for each Mb a uniformly distributed random number is chosen between 0 and 1 and the Mb is moving if > ! . Hence the number !"# controls the standard deviation of Mb mean velocities. The value !"# has been calibrated about 0.3 but it can depend on the simulation when the aim is to simulate an experiment.

Trajectory analysis tools
In this section we analyze a set of trajectories by writing them as complex linear combinations of a few basic curves. That is mathematically equivalent to performing a Principal Component Analysis (PCA) on the set of Mb displacement steps, but here it is more intuitive. Keeping in mind that a plane trajectory defined by N points defines a complex vector of dimension N, a set of trajectories is in fact a set of complex vectors in dimension N. The basic trajectories are built using a complex "Singular Value Decomposition" of this set of vectors, a mathematical tool close to a PCA and which can be outlined as follows: the first basic trajectory is the trajectory which, after a similitude (i.e. a rotation and a dilatation), best fits all the trajectories, in the mean square sense. After subtracting this first approximation from all trajectories, we obtain a second set of trajectories. The second basic curve is built in the same way as the first one from this second set of trajectories, and so on. In Figure 4, the application of this analysis for the trajectory of the experiment WT2: the first basic curve is close to a straight line and shows the direction of the trajectory (and its start-end length). The second basic trajectory shows a simple oscillation of the trajectory around this line (that can be also away and return). The third basic trajectory shows a double oscillation and so on (the basic trajectories of Figure 4A have been scaled by a weight proportional to their average contribution to fit the trajectories). This analysis is similar to a Fourier series decomposition but with basic functions which are the best adapted to each set of trajectories. Actually, less than 7 basic trajectories are needed to represent all trajectories with an error less than 2%. To sum up, for each set of trajectories an optimal set of a few basic trajectories is calculated such that each trajectory is a complex linear combination of these basic trajectories. Since a multiplication by a complex number is a similitude (i.e. a rotation and a dilatation), meaning that each trajectory is a linear combination of trajectories similar to the basic ones. Figure 4 shows the rebuilding of a set of experimental trajectories taking account 1, 2 or 3 basic functions. Now, to compare two sets of trajectories by the above analysis, basic trajectories are extracted from one set and both sets of trajectories (experimental and simulated) are expanded into this basis. Note that in order to make a meaningful comparison between two sets of trajectories, some preprocessing must be applied to the data. First the time interval must be the same for the two sets; hence we reduce the longest to the smallest. Then, since a trajectory is in fact a set of points, i.e. a complex vector, all the trajectories must have the same number of points. That is why all the trajectories are interpolated with 100 points to prevent the loss of the precision.

Crossing number assessment
The aim is to assess the number of crossings of two Mb trajectories, or to be more precise, the number of intersections of the two curves defined by the trajectories. Two irregular curves can have many intersections without a clear meaning of their number of intersections due to the lack of precision of experimental data and calculations. That is why we consider it more meaningful to count the number of pairs of curves that have any intersections. In short we do not count all crossings between two trajectories but only those occurring during the memory effect times and we count at most one crossing between two trajectories in order to avoid that two neighbouring trajectories lead to a lot of crossings. The following method is used: the domain of study (a square) is divided into small squares and we count the number of pairs of curves that are going through a small square, but only once. Notice that this method does not distinguish genuine crossings from one-sided contacts.

Justification of method choices and parameter assessment
The geometric parameters (quasi random polygon centers, the ratio of randomness, and the minimal edge length ) are chosen to mimic the biological data, here mainly the Kc polygonal shape (see Figure 2A). The number of Mb and the density of Kc are the same as in the experiments (about 220 Mb and 4,000 Kc with an average edge size of 10µm).
The main choices are in the Kc potential definition: -The maximum !"# of the Kc potential depends of the Kc polygon size but after some tests we choose to leave it invariable and since its value has no meaningful effect on the trajectories, it is fixed to !"# = 1. -In order to get a good match between experiment and simulation, the Kc attracting potential must combine suitable spatial and time variations: both are obtained with a quick potential decreasing for bound Kc and a suitable diffusion for the potential. With other choices the potential is either too flat, leading to spatial concentrations of Mb, or too irregular, leading to too random-like trajectories. The speed of increasing ! of a free Kc potential, the speed of decreasing ! of a bound Kc potential and the speed of diffusion of the factors associated with these potentials are strongly correlated with these qualitative properties and after some qualitative and quantitative tests, we calibrated the values ! = 0.2, ! = 10 and ! = 2 and the coefficient weight for the diffusion ! = 1 in order to get a noticeable and spatially consistent potential variation after a few time steps. -Note that the attractive and repulsive potentials ! and ! are of the same order of magnitude, hence the ratio ! between the Kc potential and the Mb repulsion potential is of the order of magnitude 1. Smaller or larger values lead to neglecting the effect of one of the potentials. The results are not very sensitive to the exact value. The potential ! favouring a direction has a small coefficient, between 0 and 0.15, depending of the experiment to match.       Table S1. Values and parameters used in the mathematical model