A dynamic, spatially periodic, micro‐pattern of HES5 underlies neurogenesis in the mouse spinal cord

Abstract Ultradian oscillations of HES Transcription Factors (TFs) at the single‐cell level enable cell state transitions. However, the tissue‐level organisation of HES5 dynamics in neurogenesis is unknown. Here, we analyse the expression of HES5 ex vivo in the developing mouse ventral spinal cord and identify microclusters of 4–6 cells with positively correlated HES5 level and ultradian dynamics. These microclusters are spatially periodic along the dorsoventral axis and temporally dynamic, alternating between high and low expression with a supra‐ultradian persistence time. We show that Notch signalling is required for temporal dynamics but not the spatial periodicity of HES5. Few Neurogenin 2 cells are observed per cluster, irrespective of high or low state, suggesting that the microcluster organisation of HES5 enables the stable selection of differentiating cells. Computational modelling predicts that different cell coupling strengths underlie the HES5 spatial patterns and rate of differentiation, which is consistent with comparison between the motoneuron and interneuron progenitor domains. Our work shows a previously unrecognised spatiotemporal organisation of neurogenesis, emergent at the tissue level from the synthesis of single‐cell dynamics.

Multi-parametric exploration of HES5 inter-cellular and intra-cellular dynamics used to investigate microcluster occurrence (a) and spatial periodicity (b). Panels represent maps of intercellular time delay (x-axes), repression threshold (y-axes), and Hill coefficient (plotted as separate maps in each column). Panels arrayed horizontally correspond to distinct HES5 autorepression parameter sets marked as C15, C10 and C5 with 3 repeats (e.g. C15-1 to 3) corresponding to low coherence (15% to 5%) dynamics where low coherence indicates reduced oscillatory activity in the population. Parameter sets used are presented in Appendix.       Appendix Table S1. Quantification of number of kymographs per independent experiment showing spatial periodicity following an auto-correlation test with >95% statistical significance.
Quantification of drift in Y-axis (Dorsal to Ventral) from single cell data corresponding to kymograph experiments. Quantification of displacement in apico-basal (AB) and dorso-ventral (DV) axis from single cell tracking data reporting sample size per experiment, mean of single-cell root mean squared (RMS) distance and interquartile range (IQR) over 2.5h.

Condition
Average displacement in Y-axis per 12h  Figure S6.     Table S5. smiFISH probes weaker amplification product for the second reaction which we could not re-amplify with high fidelity polymerases. Pup #3 was bred forward with a WT C57BL6/J mouse, germline transmission confirmed through PCR and sequencing, and a colony established.

Stochastic model of auto-repression coupled between cells
We adapted a HES5-parameterised single-cell model from previous work into a Notch-Deltacoupled multicellular model (Figure 5a) (Manning et al., 2019). The single cell network is the basic unit of the model and consists of negative feedback loop of HES5 protein onto its own mRNA expression described by stochastic delay differential equations (SDDEs) capable of producing stochastic autonomous oscillations. The SDDEs implement a Langevin approach over a Gillespie algorithm in favour of reduced computational cost and fixed time steps which enable easier addition of time delays (Gillespie, 2000).
To extend the single cell model to a multicellular one, the cells were coupled together with an interaction representative of the Notch-Delta pathway and its interaction with HES5. Instead of modelling every reaction in the chain of events that constitute the Notch-Delta pathway, we used a highly adjustable coupling function -an inhibitory Hill function -that requires only two parameters to specify its shape. This drastically reduced the potential number of parameters required to describe the Notch-Delta pathway that couples HES5 dynamics between cells while maintaining the ability to describe a wide variety of possible coupling function shapes that may represent the biological system. The Hill function in this model therefore describes how the production rate of HES5 in one cell is inhibited by high levels of HES5 expression in a neighbouring cell and has the form of where ( )is the mRNA production rate response of a cell to the concentration of protein in a neighbouring cell, ! is the repression threshold (also known as the effective concentration EC50) and is the Hill coefficient. ! and define the shape of the Hill function (illustrated in Figure 5b) and thus how a cell will respond to the levels of HES5 in its neighbouring cells.
Specifically, we refer to repression threshold as the inverse of the coupling strength as this parameter defines at what concentrations of HES5 in a neighbouring cell a receiving cell becomes repressed. The Hill coefficient is also capable of changing the coupling strength but to a lesser extent as it affects the slope of the curve around the value of the repression threshold, creating a more step-like response at high Hill coefficient value and reproducing a Michaelis-Menten curve at a value of 1.
For the analysis in this paper, we chose neighbours to be taken as the closest or first rank neighbours, which in the case of hexagonal geometry implies a cell can be coupled to the dynamics of a maximum of 6 immediate neighbours (illustrated in Figure 5c). The incoming protein concentration that a cell experiences is defined as the average of its neighbours where is the number of neighbouring cells and " is the protein concentration in neighbouring cell . Equation (2) assumes that a cell receives equal contribution from each of its neighbouring cells, and there is no efficiency of signalling or scaling parameter in (2) (2)) rather than a static multiplication term reduces the effects of different patterning at the boundaries of the simulated tissue and was found to result in similar dynamics as when periodic boundaries were used.
Finally, to capture the dynamics of the Notch-Delta pathway more accurately, a time delay was added to the Hill function. The use of a nonlinear function to describe the multiple steps neglects the time-delays associated with the reaction, transport/diffusion, and synthesis, so an explicit time delay parameter was added. The full form of the coupling function taking into account the neighbours and time delay is given by where the subscript indicates intercellular parameters associated with Notch-Delta signalling and $% is the intercellular time delay.
The full set of equations that describe the dynamics in an individual cell within the multicellular model is the single cell model from previous work (Manning et al., 2019) but with the introduction of (3) as a multiplicative term that affects the production rate of mRNA Where is HES5 mRNA concentration, !,% are the degradation rate of HES5 mRNA and protein respectively, !,% are the transcription and translation rates of HES5, !,% are Gaussian white noise, and , ( − " )/ is an inhibitory Hill function that describes Hes5 autorepression where &! is the repression threshold for HES5 acting on its own promoter and & is the corresponding Hill coefficient. The multiplicative nature of the two Hill functions in (4) is based on previous literature (Chikayama, 2014;Lewis, 2003) In (4) and (5), there is a deterministic part which represents the overall increase or decrease in protein or mRNA to be expected at given concentrations and then a stochastic part that accounts for random binding or non-binding of transcription and translation machinery in the cell. The stochastic part is scaled by the square root of the number of events that can happen in a given time and so the stochastic part becomes increasingly significant at low protein/mRNA numbers. The Euler-Maruyama method was used to numerically calculate solutions to the SDDEs in (4) and (5) using a time step size of 1 minute.
Parameter values for simulations of the model included in Figure 5 and Appendix Figures S5 are given in Appendix Table S4 marked as Main. We also investigated several other parameter combinations (Appendix Table S4 C15/10/5-1 to 3 and Appendix Figure S6) previously identified to reproduce HES5 tissue data by Bayesian inference (Manning et al., 2019).

Synthetic kymograph and spatial data generation
The 15µm strips used for spatial analysis in the experimental data (Figure 2a) were found to correspond to a width of 1-2 nuclei. Therefore, to obtain similar data from the model, spatial signals were generated by averaging 2 adjacent columns of cells. This single column of averaged data was used for generating each time point in the kymographs (Figure 5g,h) and also to test cluster occurrence and spatial periodicity (Appendix Fig S5c,d and Appendix Fig   S6a,b).

Parameter space explorations and parameter selection
Parameter space 2D maps were used to visualise regions of different model behaviour ( Figure   5, Appendix Figure S5 and Appendix Figure S6). These parameter spaces have intercellular time delay plotted along x-axis and repression threshold value along the y-axis and at each point an output value of the model is plotted (e.g. temporal period). For each grid point in parameter space, the average output from 5 simulations with random initial conditions is shown. Parameter selection (Figure 5f) was used to identify suitable ranges for intercellular time delay and repression threshold by comparing statistics from synthetic data against experimental statistics obtained from tissue data whereby synthetic data values found within 2.4 standard deviations away from the experimental mean were accepted. Parameter selection was performed separately for temporal period and Kuramoto order parameter (see Phase reconstruction). Temporal period was computed from synthetic timeseries of HES5 as inverse of dominant frequency peak from power spectrum reconstruction by FFT. Temporal period used for parameter selection was average period per simulation where single cell period estimates above 10h were excluded as non-oscillatory.

Rate of differentiation from synthetic data
We explored how absolute sensing of HES5 by downstream targets might affect the rate of differentiation at different intercellular repression thresholds (coupling strengths) by assuming a simple linear differentiation condition where lower HES5 levels increase the probability of a cell committing to eventual differentiation (Figure 6a). This assumption is based on the fact that HES5 represses downstream pro-neural targets which promote a cell to head towards a differentiated state.
Because the level of expression changes as the coupling strength is varied, no set absolute threshold is assumed, instead a differentiation threshold is determined from the mean population expression in each simulation. Above this threshold, a cell has a zero chance of differentiation, whereas below the threshold has a linearly increasing probability of differentiation the further it drops below the threshold. A probability is calculated at each time point with no memory of past HES5 dynamics and is calculated as where is the differentiation rate (with units of '# as probability of differentiation is calculated on every time step), and ()*+,) is the differentiation threshold.
Cells are simply marked as having differentiated in this setup, and the dynamics of an individual cell remain the same after a differentiation event, which can be thought of an asymmetric division, where the differentiating cell leaves the population and it's space in the grid is filled with a non-differentiating cell. Differentiation rate is calculated by taking the total number of differentiation events in a simulation and dividing by the total time. Single cell parameters used are shown in Appendix Table S4, and the multicellular parameters used were $% = 4, τ $% = 150 minutes, and range of repression thresholds were explored 500 < $%! < 30000.

Synthetic microcluster occurrence and diameter calculation
To detect presence of clusters of similar expressing cells in the synthetic data, spatial kymograph data was first smoothed using a polynomial fit. At each time point, the maximal peak in the polynomial was identified, followed by identification of the nearest trough below 80% the value of the peak which is consistent with a 1.2x spatial fold change in amplitude.
The presence of a peak and trough satisfying these conditions is indicative of a cluster at a specific timepoint. The frequency of cluster occurrence was calculated as the number of time points in a simulation a cluster was identified.
The microcluster diameter was calculated from spatial data by searching for groups of cells that exhibit local in-phase characteristics. We performed phase reconstruction followed by