%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ELIFE ARTICLE TEMPLATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PREAMBLE 
\documentclass[9pt,lineno]{elife}	% \documentclass[twocolumn,showpacs, superscriptaddress]{revtex4-1}
% Use the onehalfspacing option for 1.5 line spacing
% Use the doublespacing option for 2.0 line spacing
% Please note that these options may affect formatting.

\usepackage{lipsum}			% Required to insert dummy text
\usepackage[version=4]{mhchem}
\usepackage{siunitx}
% \DeclareSIUnit\Molar{M}
\usepackage{graphicx}
\usepackage{placeins}		%floatbarrier
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{wasysym}
\usepackage{dsfont}
\usepackage{color,soul}
\usepackage{dcolumn}		% Align table columns on decimal point
\usepackage{bm}			% bold math

\allowdisplaybreaks			% allow page break equations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ARTICLE SETUP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Segmentation of the zebrafish axial skeleton relies on notochord sheath cells and not on the segmentation clock}

\author[1,2,3]{L. Lleras-Forero}
\author[4]{R. Narayanan}
\author[5,6,7]{Luis G. Morelli}
\author[3\authfn{3}]{L. F. A. Huitema}
\author[1\authfn{4}]{M. VanBergen}
\author[3]{Alexander Apschner}
\author[3]{J. Peterson-Maduro}
\author[3]{I. Logister}
\author[8]{G. Valentin}
\author[4,8,9\authfn{1}]{A. C. Oates}
\author[1,2,3\authfn{1}*]{S. Schulte-Merker}

\affil[1]{Institute for Cardiovascular Organogenesis and Regeneration, Faculty of Medicine, WWU M\"unster, Mendelstrasse 7, 48159 M\"unster, Germany}
\affil[2]{CiM Cluster of Excellence (EXC-1003-CiM), M\"unster, Germany}
\affil[3]{Hubrecht Institute-KNAW and UMC Utrecht, Uppsalalaan 8, 3584 CT, Utrecht, The Netherlands}
\affil[4]{The Francis Crick Institute, 1 Midland Road, London, NW1 1AT UK}
\affil[5]{Instituto de Investigaci\'on en Biomedicina de Buenos Aires (IBioBA)---CONICET---Partner Institute of the Max Planck Society, Buenos Aires, Argentina}
\affil[6]{Departamento de F\'{\i}sica, FCEyN, UBA, Ciudad Universitaria, 1428 Buenos Aires, Argentina}
\affil[7]{Department of Systemic Cell Biology, Max Planck Institute for Molecular Physiology, Dortmund, Germany}
\affil[8]{Department of Cell and Developmental Biology, University College London, Gower Street, London, WC1E 6BT UK}
\affil[9]{Institute of Bioengineering, \'Ecole Polytechnique F\'ed\'erale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland}


\corr{schultes@ukmuenster.de}{SS-M}
% \corr{email2@example.com}{FS}

\contrib[\authfn{1}]{These authors contributed equally to this work}
%\contrib[\authfn{2}]{These authors also contributed equally to this work}

\presentadd[\authfn{3}]{Department, Institute, Country}
\presentadd[\authfn{4}]{Current address: Department of Laboratory medicine, laboratory of Hematology, Radboud University medical center. Geert grooteplein 8, P.O. 91016500 HB Nijmegen, The Netherlands}
% \presentadd[\authfn{5}]{eLife Sciences editorial Office, eLife Sciences, Cambridge, United Kingdom}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ARTICLE START
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

% \maketitle
%\begin{abstract}
%Segmentation of the axial skeleton in amniotes depends on the segmentation clock which patterns the paraxial mesoderm and the sclerotome. While the segmentation clock clearly operates in teleosts, the role of the sclerotome in establishing the axial skeleton is unclear. We severely disrupt zebrafish paraxial segmentation, yet observe a largely normal segmentation process of the chordacentra. We demonstrate that axial entpd5+ notochord sheath cells are responsible for chordacentrum mineralization, and serve as a marker for axial segmentation. While autonomous within the notochord sheath, entpd5 expression and centrum formation show some plasticity and can respond to myotome pattern. These observations reveal for the first time the dynamics of notochord segmentation in a teleost, and are consistent with an autonomous patterning mechanism that is influenced, but not determined by adjacent paraxial mesoderm. This behavior is not consistent with a clock-type mechanism in the notochord.
%\end{abstract}

\section{Theory}
\label{sec:theory}
%
The theory describes notochord sheath cells pattern formation in terms of a one dimensional reaction diffusion system with two components, an activator U and an inhibitor V, see Figure 7 - figure supplement 9A.
%
The concentration $U$ = [U] of the activator is reflected in Entpd5 concentration. 
%
The concentrations of both the activator and inhibitor species $U(x,t)$ and $V(x,t)$ depend on position $x$ and time $t$.
%
We propose a variant of the FitzHugh-Nagumo (FHN) model (32) % reference to Murray book 
%
\begin{eqnarray}
	\frac{\partial U}{\partial t}  &=& D_U \frac{\partial^2 U}{\partial x^2} + k_1 U - k_3 U^3 - k_4 V + k_0	\label{eq.pdeu}	\\
	\frac{\partial V}{\partial t}  &=& D_V \frac{\partial^2 V}{\partial x^2} + k_5 U - k_6 V				\label{eq.pdev}
\end{eqnarray}
%
where $D_U$ and $D_V$ are diffusion coefficients for $U$ and $V$ respectively, and $k_i$ are rate constants.
%
The choice of the FHN model is based on its simplicity. 
%
There are positive linear terms for the activator and negative linear terms 
for the inhibitor in both Eqs.~(\ref{eq.pdeu}) and (\ref{eq.pdev}), and 
a single nonlinearity, the cubic term for the activator that limits growth and allows for the stabilization of steady states.
%
The model is meant to represent a plausible mechanism rather than specifying the interactions of particular molecules.


% \subsubsection{Dimensionless theory}
% \label{sec:dimensionless}
%
We reduce the number of parameters by transforming this theory to a dimensionless form.
%
We first set the source term $k_0=0$ since as discussed below we need to reproduce a sequential pattern formation.
%
We introduce a lengthscale $L$ that we will take as the system size,
and a timescale $T$ and concentration scale $U_0$ to be set below.
%
In terms of these scales we define new variables $x'$, $t'$, $u$ and $v$ such that
%
\begin{eqnarray}
	x &=& L \, x'	\\
	t &=& T \, t'	\\
	U &=& U_0 \,u	\\
	V &=& U_0 \,v
\end{eqnarray}
%
and replace in the reaction diffusion equations above dropping the primes for notational convenience
%
\begin{eqnarray}
	\frac{U_0}{T} \frac{\partial u}{\partial t}  &=& \frac{D_U U_0}{L^2} \frac{\partial^2 u}{\partial x^2} + k_1 U_0 u - k_3 U_0^3 u^3 - k_4 U_0 v		\\
	\frac{U_0}{T} \frac{\partial v}{\partial t}  &=& \frac{D_V U_0}{L^2} \frac{\partial^2 v}{\partial x^2} + k_5 U_0 u - k_6 U_0 v .
\end{eqnarray}
%
We multiply both equations by $T/U_0$ to render them dimensionless
%
\begin{eqnarray}
	\frac{\partial u}{\partial t}  &=& \frac{D_U T}{L^2} \frac{\partial^2 u}{\partial x^2} + k_1 T u - k_3 T U_0^2 u^3 - k_4 T v		\label{eq.pdeudimless} \\
	\frac{\partial v}{\partial t}  &=& \frac{D_V T}{L^2} \frac{\partial^2 v}{\partial x^2} + k_5 T u - k_6 T v	.				\label{eq.pdevdimless}
\end{eqnarray}
%
Multiplying the inhibitor equation by $k_1/k_5$ and rearranging terms 
%
\begin{eqnarray}
			\frac{\partial u}{\partial t}  &=& \frac{D_U}{k_1 L^2} (k_1 T)\frac{\partial^2 u}{\partial x^2} + (k_1 T) u - \frac{k_3 U_0^2}{k_1} (k_1 T) u^3 - \frac{k_4}{k_1} (k_1 T) v		\\
\frac{k_1}{k_5}	\frac{\partial v}{\partial t}  &=& \frac{D_V}{k_5 L^2} (k_1 T)\frac{\partial^2 v}{\partial x^2} + (k_1 T) u - \frac{k_6}{k_5} (k_1 T) v
\end{eqnarray}
%
where we have highlighted dimensionless groups in parentheses.
%
We now select a timescale by setting 
\begin{equation}
k_1 T \equiv 1
\end{equation}
%
and a concentration scale by setting
\begin{equation}
\frac{k_3 U_0^2}{k_1} \equiv 1
\end{equation}
%
and define dimensionless parameter groups
%
\begin{equation}
\tau \equiv \frac{k_1}{k_5}, \qquad a \equiv \frac{D_U}{k_1 L^2}, \qquad b \equiv \frac{D_V}{k_5 L^2}, \qquad 
	\kappa_4 \equiv \frac{k_4}{k_1}, \qquad \kappa_6 \equiv \frac{k_6}{k_5}.
\end{equation}
%
With these definitions
%
\begin{eqnarray}
	\frac{\partial u}{\partial t}  &=& a \frac{\partial^2 u}{\partial x^2} + u - u^3 - \kappa_4 v		\\
\tau	\frac{\partial v}{\partial t}  &=& b \frac{\partial^2 v}{\partial x^2} + u - \kappa_6 v	.
\end{eqnarray}
%
We additionally set $\kappa_4=1$ for simplicity and we call $\kappa_6 = d$
%
\begin{eqnarray}
	\frac{\partial u}{\partial t}  &=& a \frac{\partial^2 u}{\partial x^2} + u - u^3 - v	\label{eq.fhnu}	\\
\tau	\frac{\partial v}{\partial t}  &=& b \frac{\partial^2 v}{\partial x^2} + u - d v	.	\label{eq.fhnv}
\end{eqnarray}
%
In the rest of this work we consider this dimensionless form of the theory.
%
Here $a$ and $b$ are dimensionless scaled diffusion coefficients of the activator and inhibitor species respectively, 
$\tau$ is a relative timescale, and $d$ is a dimensionless degradation constant of the inhibitor. 



We first consider the homogeneous system $\partial_{xx}u = \partial_{xx}v = 0$.
%
The resulting equations for the local reactions are
\begin{eqnarray}
	\dot u &=& u - u^3 - v	\\
\tau	\dot v &=& u - d v	.
\end{eqnarray}
%
where dots denote time derivatives.
%
Introducing functions
%
\begin{eqnarray}
f(u,v)  &=& u - u^3 - v		\\
g(u,v) &=& \tau^{-1} u - \tau^{-1} d v	,
\end{eqnarray}
%
the nullclines of the system, defined by setting $f(u,v)=0$ and $g(u,v)=0$, are the curves in the $(u,v)$ plane
%
\begin{eqnarray}
v &=& u - u^3	,	\\
v &=& d^{-1} u	.
\end{eqnarray}
%
%\begin{equation}
%v = u - u^3
%\end{equation}
%and
%\begin{equation}
%v = d^{-1} u
%\end{equation}
%
The first one is an inverted cubic that goes through the origin and the second one is a linear function with slope $d^{-1}$ 
controlled by the single bifurcation parameter $d$, see Figure 7 - figure supplement 9B.
%
Intersections of these two curves are the solutions to
%
\begin{equation}
u (u^2+d^{-1}-1) = 0
\end{equation}
%
and define the fixed points of the system where $\dot u = \dot v = 0$.
%
There is always a solution $(u_0,v_0)=(0,0)$ and for $d >1$ there are two additional solutions $u_{\pm}$ satisfying
%
\begin{equation}
u_{\pm}^2 = d^{-1}-1 .
\end{equation}
%
The linear stability of fixed points is determined by the matrix 
%
\begin{equation}
A =	\left(
	\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
	f_u  & f_v \\
	g_u & g_v \\
	\end{matrix}
	\right)
\end{equation}
%
%\begin{equation}
%A =	\left(
%	\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
%	\partial_u f  & \partial_v f \\
%	\partial_u g & \partial_v g \\
%	\end{matrix}
%	\right)_{(u_*,v_*)}
%\end{equation}
%
where 
\begin{eqnarray}
f_u  &=& \partial_u f(u,v)	\\
f_v  &=& \partial_v f(u,v) 	\\
g_u &=& \partial_u g(u,v)	\\
g_v &=& \partial_v g(u,v)
\end{eqnarray}
%
and derivatives are evaluated at the fixed point $(u_*,v_*)$.
%
The condition for stability is that 
%
\begin{equation}
\mbox{det} A =  f_u g_v - f_v g_u > 0
\end{equation}
%
and
%
\begin{equation}
\mbox{tr} A =  f_u +  g_v < 0.
\end{equation}
%
For the fixed point $(u_0,v_0)=(0,0)$ we obtain
%
\begin{equation}
A =	\left(
	\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
	1  & -1 \\
	\tau^{-1} & -d \tau^{-1} \\
	\end{matrix}
	\right)
\end{equation}
%!TEX encoding = UTF-8 Unicode
with determinant and trace
%
\begin{equation}
\mbox{det} A = (1 - d) \tau^{-1}
\end{equation}
%
%
\begin{equation}
\mbox{tr} A =  1 - d \tau^{-1}.
\end{equation}
%
Given that $\tau,d>0$ this implies that the origin $(0,0)$ is a stable fixed point if 
%
\begin{equation}
d < 1  \quad \mbox{and} \quad \tau < d .
\end{equation}
%
% Thus, as $d$ increases from low values to high values, a bifurcation occurs at $d=1$.


We consider in the following a situation in which the fixed point $(0,0)$ is stable, setting the dimensionless timescale $\tau = 0.1$ and degradation $d=0.5$.
%
Under appropriate conditions, the dimensionless reaction diffusion theory Eqs.~(\ref{eq.fhnu}) and (\ref{eq.fhnv}) can give rise to pattern formation.
%
In particular we require that the activator diffuses slower than the inhibitor, $a < b$, and here set $a=10^{-3}$ and $b=10^{-2}$.
%
Because of the signs of the derivatives near the fixed point, the type of pattern predicted is as displayed in Figure 7 - figure supplement 9C,
where there is coexistence of activator and inhibitor (32).


In the presence of random perturbations distributed along the notochord the homogeneous state loses stability due to differential diffusion.
%
A pattern may form out of this initial random background fluctuation, with segments forming almost simultaneously all along the notochord. 
%
Although segment formation is robust in this scenario and can accommodate a broad range of wavelengths,
this is at odds with the experimental observation that ENTPD5 segments form sequentially from anterior to posterior.  


We conjectured that one way to obtain a sequential segment formation is to start with an initial perturbation localized at the anterior, 
and vanishing concentrations all across the rest of the notochord. 
%
Since the anterior of the vertebrate axis is always more developmentally advanced than the posterior, 
such an anterior perturbation is a plausible hypothesis. 
%
A vanishing concentration along the notochord is important to ensure that patterning is not triggered until the wave of activator and inhibitor arrives at a given point.
%
Therefore, we start simulations with initial conditions that have zero concentration for both $u$ and $v$ across the whole domain $x \in (0,L)$,
except for a small perturbation near the origin $x = 0$.
%
The form of the initial condition is a smooth step
%
\begin{eqnarray}
u(x,0) &=& \frac{u_0}{2} \left(1 - \tanh( 5.0 (x - 0.5) ) \right)	\\
v(x,0) &=& \frac{v_0}{2} \left(1 - \tanh( 5.0 (x - 0.5) ) \right)
\end{eqnarray}
%
with a steepness 5.0 and width 0.5.
%
Initial values $u_0$ and $v_0$ are determined randomly from a uniform distribution in the interval $(0.1,0.2)$, 
see examples in Figure 7 - figure supplement 1-5.
%
Starting from such a small perturbation in the anterior,
we observe that the system is able to form a pattern sequentially, from anterior to posterior, 
see Figure 7A, Figure 7 - figure supplement 1 and Movie A.
%
Thus, the theory proposed is capable of autonomous patterning of the notochord in the absence of input from the segmentation clock.


We next introduce the effect of the segmentation clock input into this otherwise autonomous patterning system 
as a spatial dependent degradation profile of the inhibitor 
%
\begin{eqnarray}
	\frac{\partial u}{\partial t}  &=& a \frac{\partial^2 u}{\partial x^2} + u - u^3 - v		\label{eq.fhnus}	\\
\tau	\frac{\partial v}{\partial t}  &=& b \frac{\partial^2 v}{\partial x^2} + u - d v - s(x) v .		\label{eq.fhnvs}
\end{eqnarray}
%
This \emph{sink profile} $s = s(x)$ for the inhibitor has peaks at given positions along the $x$ axis, 
describing the cues that the notochord patterning mechanism receives from myotomes. 
%
At positions where $s(x)$ is large, the inhibitor is locally degraded at a larger rate. 
%
%
Note that there is no source term for the activator since we set $k_0=0$ above.
% 
This feature together with the choice of sinks instead of sources to describe the segmentation clock cues are 
motivated from the observation that Entpd5 segments form sequentially. 
%
The presence of sources would render pattern formation non sequential.


The sink profile $s(x)$ is characterized by a sink strength $S_0$, a wavelength $\lambda$ and sink wavelength variability $\sigma$. 
%
The first sink is positioned at $\lambda / 2$
and the positions $X_i$ of consecutive sinks are determined by the wavelength $\lambda$ 
with an error drawn from a uniform distribution of width $\sigma$.
%
The sink profile is built from a combination of $\tanh(...)$ functions to produce smooth peaks of steepness $\alpha$ and width $\delta S$
%
\begin{equation}
s(x) = \frac{S_0}{2} \sum_{i}^{} {
\left(
	-\tanh \left( \alpha (-X_i + x - \delta S) \right) + 
	\tanh \left( \alpha (-X_i + x + \delta S) \right) 
\right)	}	.
\end{equation}
%
In this work we fix the values $\alpha = 100$ and $\delta S = 0.05$.
%
The values of $S_0$, $\lambda$ and $\sigma$ are changed to describe the different conditions, 
see examples in Figure 7 and Figure 7 - figure supplement 6.


We consider a system size $L$ that we set to $L = 17.1$ so that the wildtype condition 
makes $30$ segments with the sink profile natural wavelength $\lambda=0.57$. 
%
We normalize axes length scales to this value in all plots.
%
For simplicity we assume that the activator and inhibitor are restricted to notochord sheath cells and 
we specify Neumann boundary conditions, that is derivatives at both ends are zero
%
\begin{equation}
\left.\frac{\partial u}{\partial x}\right|_{x=0}  = \left.\frac{\partial v}{\partial x}\right|_{x=L} = 0 \, .
\end{equation}


As described above, here we ensure the sequential character of the patterning through an initial perturbation at the anterior and 
vanishing concentrations across the notochord.
%
One may query the robustness of such scenario, since noise across the notochord hampers sequential patterning.
%
An alternative hypothesis would be to postulate an additional wavefront that propagates through the notochord 
progressively turning on the reaction diffusion mechanism of Eqs.~(\ref{eq.fhnus}) and (\ref{eq.fhnvs}) as it goes.
%
% For example the sink profile might be a function of both space and time $s(x,t)$.
%
To illustrate this we consider an alternative dimensionless form of Eqs.~(\ref{eq.pdeu}) and (\ref{eq.pdev}).
%
Turning back to Eqs.~(\ref{eq.pdeudimless}) and (\ref{eq.pdevdimless}) we select a timescale and concentration scale setting
\begin{equation}
\frac{D_U T}{L^2} \equiv 1
\end{equation}
and
\begin{equation}
\frac{k_3 U_0^2}{k_1} \equiv 1  .
\end{equation}
%
Introducing dimensionless parameter groups
%
\begin{equation}
\delta \equiv \frac{D_V}{D_U}, \qquad \gamma \equiv \frac{k_1 L^2}{D_U}, \qquad \kappa_i \equiv \frac{k_i}{k_1},
\end{equation}
%
and setting for simplicity $\kappa_4 = \kappa_5 = 1$ and $\kappa_6 = \kappa$ we arrive at the dimensionless form
%
\begin{eqnarray}
\frac{\partial u}{\partial t}  &=& \frac{\partial^2 u}{\partial x^2}		+ \gamma \left( u - u^3 - v \right)			, \\
\frac{\partial v}{\partial t}  &=& \delta \frac{\partial^2 v}{\partial x^2}	+ \gamma \left( u - \kappa v \right)			.
\end{eqnarray}
%
In this alternative dimensionless form it is straightforward to decouple the reactions from diffusion by tuning the value of $\gamma$.
%
Thus, we can introduce a wavefront $\gamma(x,t)$ that moves from anterior to posterior turning on the reactions in its wake.
%
Such a wavefront could have a biological origin in a molecular maturation gradient invading the notochord from the anterior.
%
Due to very slow dynamics before wavefront arrival, this would render the patterning mechanism more robust to noise across the notochord.
%
Yet a different possibility is a scenario of patterning in a growing domain (49), % reference to Crampin 1999
%
although here the tissue where the pattern forms exists previous to the establishment of the pattern.
%
These alternatives are certainly interesting and the underlying mechanism patterning the notochord will be subject of future work.
%
However, in the absence of experimental data supporting other hypotheses, 
here we settle on the maybe more parsimonious choice of vanishing initial conditions across the notochord.


In this work we solve the partial differential equations described above using a custom python code, see Supplemental Information. 
%
We discretize space with a discretization length $\Delta x=0.01$, and time discretization is chosen as $\Delta t = 0.9 \Delta x^2 / 2$. 
%
We integrate the partial differential equations until the pattern reaches a steady state.


\clearpage
\bibliography{library.bib}

\end{document}

