Neuronal differentiation influences progenitor arrangement in the vertebrate neuroepithelium

ABSTRACT Cell division, movement and differentiation contribute to pattern formation in developing tissues. This is the case in the vertebrate neural tube, in which neurons differentiate in a characteristic pattern from a highly dynamic proliferating pseudostratified epithelium. To investigate how progenitor proliferation and differentiation affect cell arrangement and growth of the neural tube, we used experimental measurements to develop a mechanical model of the apical surface of the neuroepithelium that incorporates the effect of interkinetic nuclear movement and spatially varying rates of neuronal differentiation. Simulations predict that tissue growth and the shape of lineage-related clones of cells differ with the rate of differentiation. Growth is isotropic in regions of high differentiation, but dorsoventrally biased in regions of low differentiation. This is consistent with experimental observations. The absence of directional signalling in the simulations indicates that global mechanical constraints are sufficient to explain the observed differences in anisotropy. This provides insight into how the tissue growth rate affects cell dynamics and growth anisotropy and opens up possibilities to study the coupling between mechanics, pattern formation and growth in the neural tube.

. Determining drag coefficients that generate anisotropic growth in simulations A) Heatmap indicating the difference between simulation and experimental tissues for different values of the drag coefficients µ and µ . The absolute value of the difference in the change in DV length of the tissue over 48h plus the absolute value of the difference in the final tissue aspect ratio (AP/DV) between simulations and experimental data was taken and then this quantity was averaged over 10 simulations. Colour map shows the average difference (log scale) from 10 different simulations per point with the experimental data from E11.5 embryos. Other parameters are indicated in Table  2. Roman numerals indicate mechanical parameters used (Λ,Γ) in Fig. 3A. B) The AP/DV aspect ratio of simulated tissues at the end of the simulation (72h) for a range of µ /µ ratios. µ = 1 is held constant and µ is varied from 0.0025 to 1. The proliferation rate is 0.05h −1 and there is no differentiation. Parameter regime V is used.   . Effect of proliferation and differentiation rate on the size and shape of simulated neuroepithelial tissue A) Number of cells in simulated tissues as a function of time for different proliferation rates and no differentiation. The proliferation rates from top to bottom: 0.05h −1 (red), 0.04h −1 (blue), 0.03h −1 (green), 0.02h −1 (purple), 0.01h −1 (grey), 0h −1 (black). In A-I: the simulation starts with an initial period of (t < 26h) during which the proliferation rate is 0.05h −1 and there is no differentiation. Afterwards proliferation or differentiation is altered as indicated. In A-H parameter regime V is used. B) The number of cells as a function of time for different differentiation rates and fixed p roliferation r ate = 0 .05h −1 . T he d ifferentiation r ate i s 0 h −1 ( red), 0 .01h −1 ( blue), 0 .02h −1 ( green), 0 .035h −1 ( purple), 0 .05h −1 ( grey), 0.075h −1 (black). The parameters were chosen so that the net growth rate in B is similar to that in A for curves of the same colour (except black). C) Overlaid data from A and B. D) Aspect ratio of the tissue AP to DV extension for varied proliferation and no differentiation (same as main . Influence of the ratio of the AP/DV drag coefficients on the shape of tissue with different rates of differentiation Increasing the differentiation rate increases the AP/DV aspect ratio of simulated tissues for all µ /µ ratios larger than 1. The AP/DV aspect ratio at the end of the simulation (72h) is reported. The simulation starts with an initial period of (t < 26h) during which the proliferation rate is 0.05h −1 and there is no differentiation. Afterwards differentiation is altered: 0h −1 (purple, same as in Fig. S.3B), 0.05h −1 (red), 0.1h −1 (solid light blue). The indicated µ /µ ratio is fixed from the start. The parameter regime V is used. In the case µ /µ =1, the final AP/DV aspect ratio is <1. This is because the tissue is initialized with hexagonal mesh of 10 by 10 cells (honeycomb pattern) that has the AP/DV ratio = √3/2 ≈ 0.87. If the tissue is initialized with hexagonal mesh of 10 cells in DV direction and 12 cells in AP direction (initial AP/DV ratio ≈ 1.04), then the end result is also isotropic as expected (empty light blue dot). The reported errors are SE.   The domain of the vertex model is a cylinder, the radius, R, and height, H, of which can change in time. A natural set of coordinates for the vertex i within the cylinder are (θ i ,z i ), where θ i are cylindrical polar angles andz i = z i /H, where z i is the coordinate in the direction along the length of the cylinder, see Fig. S9. Note that in simulations we actually impose periodic boundary conditions also in the direction along the length of the cylinder. The movement of vertices is determined by an energy function, E.
We assume overdamped movement, with the drag coefficient for vertices within the cylinder (whilst the cylinder size is fixed) being µ. Therefore: In addition the cylinder itself can grow. We assume that cells (whose average areas are constant) experience some drag proportional to their speed, both radial and tangential (frictional) as the outside surface of the cylinder expands. Radially the speed of each cell is dR dt and along the length of the cylinder, the velocities are on average 1 2 dH dt relative to a fixed point at one end of the cylinder 1 . Let us define the relevant drag coefficients to be µ radially and 2µ along the length of the cylinder. The equation for the radius and the height are determined by the balance between the radial component (the height component) of the potential forces and the accumulated drag on all of the cells: where N c is the number of cells. We note that the three drag coefficients µ, µ and 2µ are likely to be different. We also note when a cell is added its energy is added to E.
The coordinates that we actually use in the vertex model code are given by (x i , z i ), where x i = Rθ i and z i = Hz i , which are effectively rectangular coordinates in a folded out cylinder. The equation for the radius, thus becomes: The magnitude of the overall force is the same if we consider the middle of the cylinder as being stationary. Therefore Similarly, In the code, we update x i and z i first according to equation (SM.1) and (SM.2), (i.e. in time step ∆t we add R∆θ i and H∆z i ), then we multiply by a factor R+∆R R and H+∆H H respectively, where we call ∆R R and ∆H H "expansion" and these are given by − ∆t

II. BIOPHYSICAL PARAMETER SPACE
We study the effect of the biophysical parameters on the cell dynamics. The energy (Eq. 4) is determined by the vertex positions and parameters K, Λ, and Γ. Magno et al. [2015] find parameter regions that yield different behaviours depending on the mechanical parameters.
The energy E for an individual cell is given by: The interfacial tension γ and pressure Π are given by A reduction in cell perimeter reduces the energy and is therefore favoured. This favours cells that take a regular shape, hence we require the interfacial tension to be positive, This should be true for the equilibrium cell size. Suppose a cell has a fixed shape and denote the square root of its area by l and its shape index by s = L/ √ A where L and A are the perimeter and the area cell. The equilibrium size is determined by setting ∂E ∂l = 0: (SM.11) thus at equilibrium γs = 2Πl, and therefore 2Kl 3 = − Λ 2 s + (2KA 0 − Γs 2 )l. (SM.12) Development: doi:10.1242/dev.176297: Supplementary information We note that when l = 0, ∂E ∂l = Λs, so the state l = 0 is stable if and only if Λ > 0. Eq. SM.12 has a single positive solution corresponding to an energy minimum if Λ < 0 (meaning cells converge to a single equilibrium area). If Λ > 0, Eq. SM.12 can either have zero positive solutions (in which case all cells collapse to zero area) or two positive solutions. In the latter case, the larger positive solution will correspond to an energy minimum and the small positive solution to an energy maximum. Thus cells with area greater than the unstable value will expand their area towards a fixed positive value, while the area of small cells will collapse to zero.
We note that the condition for the state l = 0 to be unstable is independent of the target area. In our model, the target area varies slowly compared to the movement of vertices and thus we expect this property to hold for the same parameter values in the model with varying target area too. The condition for regularly shaped cells does depend on the target area, since when interfacial tension is zero, the minimum energy will be attained when cell area is equal to target area. With other parameters fixed, as the target area increases, the interfacial tension will also increase and so cells will tend to be more regular. Thus it is possible that small cells will be irregularly shaped, but larger ones will be regularly shaped. This is because for small cells the adhesive force will dominate the force from the actomyosin ring, whereas for larger cells, the opposite will be true.
The final condition that we need to determine is the condition for there to be two positive steady states instead of zero when Λ > 0. There are two positive states if and only if This condition depends on our variable target area. What really matters, however, is whether cells are likely to grow at any size that they attain (i.e. it does not matter if we are in the region with two stable areas if the cells are so small when they are born that they shrink to zero). A necessary condition for cells to shrink from their size at birth is that ∂E ∂l l= √ Ac/2 < 0 with the target area given by its value at the start of the cell cycle. A sufficient condition is that ∂E ∂l l= √ Ac/2 < 0 with the target area given by its minimal value. Now (SM.14) The cells will thus grow if 2K A c /2(A 0 − A c /2) > Λ 2 s + Γs 2 A c /2. (SM.15) Let us use the notation of Farhadifar et al. [2007], with the A 0 that they use to create the nondimensional parameters given by its mean value, A 0 during the cell cycle in our model. The nondimensional parameters arē Λ = Λ K A 0 3/2 andΓ = Γ K A 0 . Let A c /(2 A 0 ) = c and the value of A 0 in Eq. SM.15 be β A 0 . Then cells will grow ifΛ If the model of target area variation with the cell cycle and the critical area for cell division are fixed, then, for a given shape index (e.g. that of a regular hexagon), this forms a diagonal region in the phase space.
We note that equations SM.13 and SM.15 imply that cells with lower shape indices (i.e. more regular and with larger number of sides) are more likely to be stable.
In Fig. S10, we show the phase diagram in the space of the parametersΛ andΓ. Region I is where cells are expected to have negative interfacial tensions and so have irregular shapes since they do not want to minimise their perimeters. The boundary of this domain is at Λ + Γs √ A 0 = 0, soΛ +Γs = 0. We plot this for hexagonal cells (s = 2 × 3 1/4 2 1/2 ). In Region II, cells take regular shapes and the only stable size is positive. The other boundary of this domain is at Λ = 0 or Λ = 0. In Region III, cells take regular shapes and small cells collapse, whilst larger ones grow to an equilibrium size. In Region IV, all cells collapse. The boundary between these regions is given by A 0 = Γs 2 2K + 3 4 Λs K 2/3 or 2 =Γs 2 + 3 2 (sΛ) 2/3 . We plot this for hexagonal cells. We also show as the dotted line, a line below which all cells should grow (Eq. SM.16). Here β = 0.68/1.25 is the minimum value of A 0 divided by its mean value. The minimal value is calculated from Eq. (1), assuming mean growth rate, and c = 1.3/(2 · 1.25) is the critical area divided by twice the mean value of A 0 (all nondimensional units).

III. TISSUE ASPECT RATIO DEPENDS ON EXPANSION RATE, THE VALUE OF Γ AND ON HOW EFFICIENTLY CELLS REARRANGE
To gain insight into the anisotropic growth of the tissue, we consider a slightly simpler system. We assume all cells are identical with area A and we assume that Λ/Γ << 1, so that the target perimeter is negligible. We assume cell division is unorientated. Let the dimensions of the neural tube be R(t) and H(t) and the number of cells be N c (t). Therefore N c A = RH. For ease of analysis, we assume that cells are rectangular. In the direction of R, there are n 1 cells of length l 1 , with R = n 1 l 1 and similarly in the other direction there are n 2 cells of length l 2 . Therefore, H/R = (n 2 l 2 )/(n 1 l 1 ) and N c = n 1 n 2 . We start with an equal number of square cells in each direction. We consider two extreme cases. In the first case, intercalation can take place rapidly, allowing the number of cells in each direction to change. In this case, any change in the tissue shape and the dissipation of drag forces will occur predominantly by changing the number of cells in each direction, rather than by altering cell shape. In such a scenario, the shape of individual cells will be relatively unaffected by anisotropic forces and will be square, hence n 2 /n 1 = H/R. In the second case, intercalation does not take place, in which case as the aspect ratio of the whole tissue changes, the cells change shape accordingly. In this case, the relative number of cells in each direction does not change, since cell division is unorientated and intercalation impossible, hence n 1 = n 2 and l 2 /l 1 = H/R.
The expansion of the tissue in dimensions R and H and perimeter L = 2(l 1 + l 2 ), is given by SM.3. On the shortest timescales, we assume only l 1 and l 2 change. Therefore Similarly, In the first case, µ /µ dR/dH = H/R and the change in R 2 will be µ /µ times the change in H 2 , so the aspect ratio, H/R, will tend to µ µ . In the second case,

(SM.19)
We expect that in the long term the growth will equilibrate so that the average areas and aspect ratios of cells will tend to constants, which we call r = H/R and A = RH/N c . Then we have that the length and height of the tissue satisfy If we assume that cell numbers grow exponentially with rate λ, so that dNc dt = λN c , then substituting into Eqs. SM.19 gives Dividing by √ N c A yields simultaneous equations for r and A 0 − A with solution r = µ λ + 8Γ µ λ + 8Γ A 0 − A = 8Γ + (µ λ + 8Γ)(µ λ + 8Γ) 2K .
This means that when growth is slow so that µ λ, µ λ << 8Γ, as t → ∞, the aspect ratio should tend to one. When growth is fast so that µ λ, µ λ >> 8Γ, the aspect ratio should tend to µ µ . As the exponential growth rate of the tissue varies from small to large, the asymptotic aspect ratio will vary monotonically between these two values. Slow tissue growth is defined by the regime in which perimeter forces greatly exceed drag forces and rapid tissue growth by the regime in which drag forces greatly exceed perimeter forces. Cell pressure due to area elasticity is always important. In slow growth, pressure generated by the difference between target area and mean area is balanced by the stretched actomyosin ring, whereas in rapid growth it is balanced by drag.
The most extreme aspect ratios are achieved when the actomyosin ring exerts a negligible force or when tissue rearrangement is extremely efficient.