Constraints on somite formation in developing embryos

Segment formation in vertebrate embryos is a stunning example of biological self-organization. Here, we present an idealized framework, in which we treat the presomitic mesoderm (PSM) as a one-dimensional line of oscillators. We use the framework to derive constraints that connect the size of somites, and the timing of their formation, to the growth of the PSM and the gradient of the somitogenesis clock period across the PSM. Our analysis recapitulates the observations made recently in ex vivo cultures of mouse PSM cells, and makes predictions for how perturbations, such as increased Wnt levels, would alter somite widths. Finally, our analysis makes testable predictions for the shape of the phase profile and somite widths at different stages of PSM growth. In particular, we show that the phase profile is robustly concave when the PSM length is steady and slightly convex in an important special case when it is decreasing exponentially. In both cases, the phase profile scales with the PSM length; in the latter case, it scales dynamically. This has important consequences for the velocity of the waves that traverse the PSM and trigger somite formation, as well as the effect of errors in phase measurement on somite widths.


Phase profile in a PSM consisting of a finite number of cells
In this section, we will calculate phase profiles for a system, in which a finite number of oscillators are placed on line. We assume that new oscillators are added at the left end of the line with time intervals T g = m L /T 0 , with m L ∈ N, T 0 > 0, and that the rightmost m R ∈ N oscillators are removed with time intervals T s = T 0 (it is straight forward to substitute another value of T s in our analysis). Here, T 0 is the time period of the oscillator at the left end of the line, corresponding to the posterior end of the PSM. For simplicity, we furthermore consider the case m L = m R = m, corresponding to a line in steady state: In a time interval T 0 , it grows as much as it is shortened. The system is illustrated in Fig. 1 in the main text. If we choose t = 0 to coincide with the addition of an oscillator at the left end of the line, and removal of m oscillators at the rightmost part of the line, and assume that the line consists of N oscillators at t = 0, the number of oscillators on the line is given by Because the line of oscillators changes its length with time, each oscillator will also change its position on the line, relative to the left end (the time dependent position of an oscillator starting on the leftmost position on the line is shown on the leftmost plot, Fig. 1 A in the main text. Because new oscillators are added on the left end of the line, an oscillator effectively moves one position to the right each time a new oscillator is added. Likewise, every time m oscillators are removed from the right end of the line, each oscillator, remaining on the line, effectively moves to the right relative to the length of the line. To formulate an expression for the relative position of an oscillator on the line, let us first consider the case where no oscillators are removed on the right hand side. The line only grows, and it does this by the addition of oscillators on the leftmost end of the line. In this case, no matter the length of the line, the number of oscillators to the right of oscillator i is constant. Assuming that the line initially had length L(0) = N , and that the initial position of the oscillator was 0 ≤ i 0 ≤ N − 1, we can exploit this fixed distance to the rightmost end of the line in writing down an expression for the relative position of the oscillator at time t, x only growth (t) = L(t) − (N − i 0 ) L(t) − 1 . (2) In the above expression, the numerator is an integer, which is N −i 0 smaller than the total number of oscillator on the line. Dividing by the total length (minus 1) yields a number in the interval [0, 1], the relative position of the oscillator. When including periodic removal of oscillator from the right, the expression for the relative position of an oscillator becomes This expression is similar to Eq. (2), except for the term that is added to i 0 in the numerator. This term, t T0 m accounts for the removal of m oscillators in time intervals of length T 0 . This is taken into account because every time oscillators are removed to the right of oscillator i, its distance to the rightmost end of the line changes, and this is what we use to calculate x i (t).
With the above definitions, we are almost ready to calculate steady state phase profiles. First, however, we must define how oscillation periods change as a function of position on the line of oscillators.
We will assume that oscillation period increases linearly from T (x = 0) = T 0 at the posterior end (leftmost end of the line, where new oscillators are added) to T (x = 1) = (1 + λ)T 0 at the anterior end (rightmost end of the line, where oscillators are removed from). This assumption is based on experimental observations in mouse mPSMs, as discussed in the main text. λ has been determined to be between 0.25 and 0.30 experimentally. Furthermore, we assume that newly added oscillators have initial phase identical to that of the oscillator which occupied the leftmost point on the line until the moment when this new oscillator was added. We implement this assumption by defining the period of oscillators having negative spatial positions to have a period and initial phase identical to that of the oscillator on the position x = 0 (still referred to as the leftmost oscillator, even though oscillators can now take negative spatial positions, corresponding to oscillators that have not been added to the line yet). Thus, the period distribution is, This allows us to write down the phase of oscillator i at time t, where t intro is the time at which the oscillator is introduced at the leftmost end of the line. We can calculate the steady state phase distribution of a line of minimum length N , given T g = m/T 0 . We can do this, by identifying N oscillators that will take positions 0, 1, . . . , N − 1 at some point, right after oscillators are removed from the rightmost end of the line. First, we notice that all oscillators that end up at the rightmost position, N − 1 immediately after removal of oscillators, must previously have taken positions N − 1 − mn, n ∈ N, at corresponding times. One oscillator, which has negative position, but satisfies this, is an oscillator Without loss of generality, we assume that the leftmost oscillator of the line has phase 2π at time t = 0, and denote the time at which an oscillator is introduced at the line, t intro = −x 0 T g = |x 0 |T g . Using Eq. (7), Here we split the last two integrals to make all upper boundaries coincide with oscillator removal. We now convert all integrals to sums over the positions the oscillator takes in each time interval, and insert t intro = |x 0 |T g . The formula for the phase then becomes where we have defined p max = m − (|x 0 | mod m), which is the position an oscillator with initial condition x 0 takes on the line, immediately after oscillators are removed for the first time following its addition to the line, and c max = N/m − t intro /T 0 , which is the number of oscillator removals that an oscillator with initial condition x 0 experiences after being added to the line before t = N/m T 0 is reached. In the last term, the first sum takes oscillator removals into account, while the second sum ensures that all m positions an oscillator takes between oscillator removals are counted.

Phase profile in a continuous PSM
In this section, we extend our methods from the previous section to lines of infinitely many oscillators. We will use this method to solve two different example problems. We consider a line of infinitely many oscillators.
Oscillators are constantly added on the left end of the line, and a fraction of the line length is removed from the rightmost end of the line in time intervals of T 0 . For simplicity, we assume that the length of the line grows linearly between removal of oscillators, and hence, between removals, the length of the oscillator line we define L(t) = L 0 (1 + βt/T 0 ), where βL 0 is the difference between the maximum and minimum lengths of the PSM. If we now assume that oscillators are removed at times t = nT 0 (this amounts to assuming T s = T 0 ; once again, it is straight forward to replace this value of T s with another), n ∈ N, and that this removal restores the length of the oscillator line to the length it had at t = 0, the length is described by As was the case in our discrete description of the line of oscillators, the position of each oscillator relative to the length of the line, effectively moves right as new oscillators are added at the left end of the line. To express the position of an oscillator as a function of time, we again exploit the fact that the distance between the rightmost point of the line, and an oscillator is constant between removals of oscillators. We write down the relative position of an oscillator as a function of time in the same way as we did in the previous section. First, if no oscillators are removed from the right end of the line, the relative position of an oscillator which had position Taking removal of oscillators into account means that x 0 → x 0 + β t/T 0 , since each point moves βL 0 right every time βL 0 is removed from the line from the rightmost end. Therefore, Positions on the line have values x ∈ [0, 1]. With the same reasoning as in the previous section, for a steady state phase profile, if an oscillator ends up in position x = 1 at time t = nT 0 , n ∈ N, it has occupied the same positions as all oscillators that ended up at x = 1 at t = n − T 0 , n − ≤ n − 1. From this follows that these oscillators were added to the line of oscillators at corresponding times between two removals of oscillators.
We can characterise such oscillators by the negative position they held at the final oscillator removal before they were added to the line. If the oscillator is added to the line at time 0 ≤ t intro ≤ T 0 , this negative position is x start = −βt intro /T 0 . If an oscillator ends up at position x = 1 at a time t = nT 0 , this negative position is That is, given the period gradient in Eq. (4), and that a newly added oscillator has the same phase as the oscillator to its immediate right, all oscillators that end up at x = 1 at a time t = nT 0 were added to the leftmost end of the line at time t intro = −x start T 0 /β, with phase φ intro = φ 0 + t intro 2π/T 0 , where φ 0 is the phase of the leftmost oscillator of the line right after a removal of oscillators at the right hand end of the line. We can write down the phase of an oscillator, starting at any position x 0 , at any time after it is added to the line, With this formula we can calculate the phase at time t, for specific initial conditions x(t = 0) = x 0 . Using this, we can get the steady state phase profile as it looks immediately after oscillator removal, if we calculate φ(t) for a set of initial condition that occupies x ∈ [0, 1] at a time t = nT 0 . An interval of initial conditions that satisfies this is The right boundary of this interval, we already concluded will end up at position x = 1. It will arrive at position x = 1 at time t = 1/β T 0 , and will eventually be removed from the position at t = ( 1/β +1)T 0 . The lower boundary of the interval is exactly L 0 to the left of this point, and hence will be at position x = 0, when the rightmost oscillator of the interval arrives at x = 1. This reasoning is similar to the one we used in the discrete case. For this reason, we calculate φ(t = 1/β T 0 ) for all initial conditions in the mentioned interval, which gives us the phase profile as a function of initial conditions, x 0 , right after oscillator removal at time t = 1/β T 0 , φ(x 0 , t = 1/β T 0 ). From this, we obtain the steady state phase profile after oscillator removal φ(x) by replacing x 0 → x − 1/β β. For our convenience, we will split the contribution to the phase of the oscillator in question into four parts: 1) The initial phase; 2) Phase acquired before the oscillator is added to the line; 3) Phase that the oscillator acquires after being added to the line, but before the first oscillator removal happens after this; 4) Phase that is acquired at later times. Assuming that all oscillators with negative initial position have initial phase 2π, the phase of an oscillator with position x(t = 0) = x 0 ≤ 0 can then be expressed, In this expression, the i th term corresponds to the i th contribution stated above. We will solve the three integrals one by one. The first integral has no explicit time dependence and the solutions is To solve the second integral, we must know T (t ). To know this, we need to know the position of the oscillator as a function of time. We can use the rescaling arguments given above, and write the position for t intro ≤ t ≤ t intro /T 0 T 0 as The last oscillator removal prior to t intro happens at t intro /T 0 T 0 , and the first oscillator removal following t intro happens at t intro /T 0 T 0 . For this reason, . Inserting this in the expression above, we get We can insert this in the expression for the period as a function of position in Eq. (4), and insert in the second integral above. We get with Having solved this integral, we now turn to the last integral in the expression of the phase at time 1/β T 0 . The final integral represents the phase that the oscillator acquires after the first oscillator removal. This time interval is an integer number of periods of the leftmost oscillator. The number of periods that the leftmost oscillator goes through before the time t = 1/β T 0 is reached, we denote c max = 1/β − t intro /T 0 . From this insight, we can write the integral as a sum c max integrals over time intervals of length T 0 . This is expressed as follows, To solve this we once again need to know the length of the line at the time at which the integrand is evaluated. In t + tintro T0 T 0 + cT 0 the final two terms are both integer powers of T 0 . Since L(t) = L(t + jT 0 ), j ∈ N, we know that L(t + tintro T0 T 0 + cT 0 ) = L(t ). To write down the position of the oscillator at time t + tintro T0 T 0 + cT 0 , we need to know the position of the oscillator at time tintro T0 T 0 + cT 0 . After each full time interval T 0 , the line is restored to length L 0 , and all oscillators have moved βL 0 to the right. For this reason, at time tintro T0 T 0 + cT 0 , an oscillator with initial condition x 0 has position x( tintro T0 T 0 + cT 0 ) = x 0 + t intro /T 0 β + cβ. We can use these pieces of information to express the position of the oscillator at time t + tintro T0 T 0 + cT 0 as follows This, we can insert into Eq. (4) to get the oscillator period at the time in question. The final integral then becomes The steady-state phase profile is obtained by adding all 4 contributions, and this is the black curved plotted in Fig. 2 in the main text.

Phase profile in a PSM that does not grow
In the previous sections we analyzed a line of oscillators that grows as much as it is shortened in one posterior period. We now turn to another important special case: A line in which there is no growth, only oscillator removal. In this case, the length of the line is conserved between oscillator removals, and an oscillator does not change its position relative to line length between oscillator removals. We will assume that 1) the phase difference between the two ends of the line is Φ before at oscillator removal (Inserting Φ before = 2π is the special case of mouse mPSMs with no growth); 2) That a phase profile is rescaled to line length but otherwise identical (mod 2π) at oscillator removal; 3) That oscillation period increases linearly along the line like above; 4) That T s = T 0 (once again, it is straight forward to substitute this assumption with another value of T s ; 5) That theφ rightmost phase is removed at oscillator removal -this means that a certain fraction of the line length 1 − x c is removed from the right end of the line at oscillator removal (that is, the oscillator with position x c ∈ (0, 1) before oscillator removal has position x = 1 after oscillator removal). In assumption 5), the value ofφ is intimately connected to the period-profile on the line; we will be using the experimentally observed valueφ = 0.21 · 2π. With these assumption we will estimate the phase profile over the line and determine x c , or equivalently the fraction of the oscillator population that is removed at each oscillator removal 1 − x c . Assumption 5) above means that an oscillator which has position xx c just before oscillator removal will have position x until next time oscillators are removed. The oscillator changes its phase δφ(x) = 2πT 0 /T (x) between these two consecutive oscillator removals. But because the leftmost oscillator acquires 2π of phase between consecutive oscillator removals, and the phase profile is in steady state, φ(x, t = nT 0 ) = 2π +φ(x, t = (n + 1)T 0 ), n ∈ N. These insights make us capable of writing down the following equation for the phase of the oscillator on position x just before oscillator removal In accordance with the assumptions above, we, without loss of generality, take the phase in the line endpoints to be φ(x = 0) = Φ before and φ(x = 1) = 0 just before oscillator removal. We now work towards an expressions that will allow us to determine the phase in infinitely many different points, and determine x c under our above assumptions. Evaluating Eq. (29) This is a recursive relation. Given a period distribution, which in this case is T (x) = T 0 (1 + λx), we can use this to determine x c such that Φ before and 0 are the phases of the endpoints of the line (these are specific to our example biological system and could be chosen differently if wanted). We can now insert δφ(x) = 2πT 0 /(1 + λx), and reduce the expression in Eq. (31) to obtain the recursive relation This is an equation in x c , which is nontrivial to solve analytically. Numerically, we evaluate the sum to e.g. n = 550 for different values of x c . For Φ = 2π, we find that x c = 0.767622 solves the equation. A bound for the error on this evaluation can be found by the following estimation where we used that x i c + 1/λ > 1, and . This yields an error on the evaluation smaller than 10 −15 .
Having estimated x c = 0.767622, we now know the fraction of oscillators that are removed periodically, 1−x c = 0.232378, and can plug x c = 0.767622 into Eq. (33) to obtain the steady state phase in any point x m c , m ∈ N. In the previous sections, we plotted phase profiles after oscillator removal, and not before oscillator removal like here. The steady state phase profile after oscillator removal is obtained by replacing x → x/x c for the evaluated points, and removing the point x = 1. This is plotted in Fig. 5A in the main text.

Changing somite width in a PSM that does not grow
In the previous section, we determined the somite width in a PSM that does not grow. This we did for a specific value of λ. In this section, we imagine perturbing the period gradient such that every oscillator has its period altered by an additive amount ξT 0 . So the new period distribution is T (x, ξ) = T 0 (1 + xλ + ξ). We assume that a somite is formed once every posterior period, and we assume that the phase width of the somite is equal to the phase difference that occurs between posterior and anterior in one posterior period. Lastly, we assume that the phase difference between anterior and posterior is 2π at the time of somite formation, and that somites form with period T s = T 0 . It is straight forward to substitute the assumed values of Φ before and T s with other values.
Eq. (29) is still valid, except that we now have an additional variable, δφ(x, ξ) is equal to the difference that occurs between posterior and anterior in one posterior period, and is given by We now proceed as we did in the previous section. Inserting φ(x n c ) recursively gives us the equation If we use φ(x 0 c ) = 0, we can calculate and with this, we can demand that the posteriormost oscillator has phase Φ before at somite formation, So for a given ξ, the corresponding x c satisfies the equation We evaluate the first 550 terms of this sum, for Φ before = 2π, and use this to find the x c that satisfies Eq. (43). The error on this evaluation is estimated as we did in the previous section, This holds if (1 + ξ)/λ > 1. This is definitely true for ξ ∈ [−0.5, 1], which is the range that we plot x c (ξ) for in Fig.5B in the main text. In the main text we refer to x c (ξ) as the physical somite width.

Physical somite width as a function of period perturbation size
The calculation in Supplementary Section 2 gave us the phase profile in a steady-state PSM whose maximal length is βL 0 longer than its minimal length L 0 . In this calculation, we have not assumed a specific Φ before . For different choices of β and λ, we get different values for the phase difference Φ before , and phase width φ. This lets us examine how perturbing the period for all cells in the PSM by an amount ξT 0 changes the physical somite width. We will do this with the following approach 1. Choose a parameter β, 2. Try different input values of λ, and find the value that corresponds to the wanted Φ before , e.g. Φ before = 2π, 3. Calculate what value of ξ, corresponds to the determined value of λ.
The first two steps in this approach are straight forward to carry out, using the analytical expression for the phase profile in a steady-state PSM in the continuum limit. We only need to figure out how to convert a chosen λ-parameter to a ξ-value. First, we note that if oscillators are removed once every posterior period (assuming T s = T 0 , other values for T s are easy to plug in to the calculations), the phase width of a somite in a PSM with posterior period T min and anterior period T max is given by, So the phase width is not determined by the difference between the posterior and anterior periods, but by the ratio between the posterior and anterior periods. We will now take advantage of our theoretical framework from Section 2 of this supplementary file. In our framework, we can choose any value for the difference in periods over the PSM, λ, we like. However, suppose that we know that only a single value λ = λ 0 := 0.21/0.79 corresponds to the period gradient in the PSMs observed in an unperturbed experiment. Suppose all other values of λ correspond to the period gradient in experiments where an additive perturbation ξT 0 affects all cells in the PSM. In the rest of this section, this is what we will assume. By assuming this, we will provide a formula for matching the chosen λ to a unique value of the perturbation size, ξ, given λ 0 . For a chosen parameter λ, the ratio between periods in the simulation is As described above, we now assume that any T min from our simulations can be written T min = T 0 (1 + ξ), and likewise for T max = T 0 (1 + λ 0 + ξ). If we demand that the observed ratio between periods f sim is equal to the ratio between these perturbed periods, we get Solving for ξ gives us This formula lets us match the λ-value, which ensures the wanted Φ before for a chosen β, with a ξ-value. We plot matching β and ξ-values for Φ before = 2π in Fig. 4B in the main text.

Phase profile shape and somite size for different growth conditions
In the previous sections, we have examined somite size, and phase profiles in the absence of growth in PSMs, and in PSMs with steady-state length. We found that while the phase profile was convex in the absence of growth, it was concave in PSMs with steady-state lengths. We also found that the physical somite size was larger in PSMs that do not grow. In this section, we will use perturbation theory to argue • For any T (x), x ∈ [0, 1], which is an increasing function of x, and PSM that has steady-state length, the corresponding steady-state phase distribution is concave.
• How the steady-state phase distribution could become convex, if the PSM does not have steady-state length.

Concave phase profiles if PSM length is in steady state
First, we examine the phase profile in a PSM in steady state. Suppose that the PSM consists of a larger number of cells. We use the continuous variable x ∈ [0, 1] to describe the position of each cell relative to the posterior end (at x = 0) and the anterior end (at x = 1). Let T (x) be the period gradient of the PSM, and let this be increasing from posterior to anterior. Suppose that two cells have initial positions x(t = 0) first := x 0,first equal to x 0 = x * , and x 0,second = x * + , where 0 ≤ x * < 1, and 0 < 1. Let us assume t = 0 to be immediately after somite formation, and let the phase difference between the two cells be δφ = φ(x) − φ(x + ) > 0. We now examine how the phase difference between these cells changes between t = 0, and just after the following somite formation at t = T s . The change in phase difference between the two cells in this time period is Now, since 1, we expand the fraction Here we used Eq. (12) to calculate ∂x(t)/∂x 0 . Inserting this expression in Eq. (51) yields, Since T (x) is increasing and positive, and since L(t) is positive and increasing between somite formation, ∆φ > 0. This means that the phase difference between the two cells increases between the two somite formations. The phase difference is the same after the somite formation at t = T 0 , and because the PSM length is in steady state, the difference in position between the two cells is still at t = T 0 . We can determine whether the phase profile is convex or concave by comparing whether the phase profile is decreasing more quickly at positions that are more posterior or more anterior . This tells us whether the phase profile is convex or concave because a decreasing, concave function has a negative second derivative, while the second derivative is positive for a decreasing, convex function (See Fig. 1S in this Supplementary file). The phase x φ(x) ∂ 2 x > 0 Figure 1S: Illustration of decreasing functions that are concave (red), and convex (blue). The concave function has a negative second derivative, while the convex function has a positive second derivative.
From this we conclude that the steady-state phase profile decreases faster as x is increased. Or equivalently: the steady-state phase profile is concave.

Convex phase profile with T (x) linear, increasing
In the previous subsection we found that the phase difference between two cells increases between somite formation if T (x) is an increasing function. This was expressed in Eq. (51). With this in mind, one may wonder how one can obtain a convex phase profile with an increasing period profile T (x), as we found in the case of no growth. The answer lies in the shortening of the PSM: Even though the phase difference between two cells increases between somite formations, the somite formation itself causes the difference in position of the cells to grow relative to the length of the PSM. This may influence the ratio between the phase profile gradients greatly. Let us return to the case of a PSM with period-profile T (x) = T 0 (1 + xλ), that does not grow. We can evaluate Eq. (51) in this case, If two cells were a distance apart from each other before somite formation, and if x c is the length of the PSM after somite formation, the distance between the same cells after the somite formation will be /x c relative to the new PSM length. Taking this into account, we can now calculate the ratio between the phase profile gradient between two cells at two consecutive somite formation events, as we did in the previous subsection, δφ +∆φ /xc δφ = x c 1 + 2πλ δφ (1 + λx 0 ) 2