Macroscopic, semi-quantitative model of a two-scale process
The macroscopic spatio-temporal model used to describe the patterning process duringthe MHB formation considers two biological scales. The single cell scale, on which asemi-quantitative model for the transcription factor activation is employed, and thetissue scale, on which morphogen concentrations and gradients are regarded. The stateof the individual cells is described by the activities of the four transcriptionfactors, Otx2, Gbx2, En and Pax and the twomorphogen encoding genes Fgf8, Wnt1, which are described in theBackground section, and the morphogens Fgf8 and Wnt1. The interactions areschematically displayed in Figure 1A. The joint dynamicsof transcription factor activitiesu(t,x) = (u1(t,x),…,u6(t,x))Tand morphogenv(t,x) = (v1(t,x),v2(t,x))Tare described by an eight-dimensional system of partial differential equations (PDEs)on the spatial domainx ∈ U = [ − L,L].As the up- and down-regulation of transcription depends only on the concentration oftranscription factors within the individual cells and the local concentration ofmorphogens, the time evolution of the transcription factor activity is governed by areaction equation,
with production term and degradation rate . While and are constant,B
i
(x,u(t,x),v(t,x)) is potentially anonlinear function of x, u(t,x) andv(t,x). B
i
represents theHill-type regulation of transcription factor i by transcription factors andmorphogens (Figure 1C). The initial conditions for thetranscription factor activity are given by the spatial functionsu0i(x),i = 1,…,6, which describe quantitatively the expressionpattern at E8.5 (Figure 2B upper panels). Hence,
(1)
We assumed that the morphogens are produced proportionally to the activation of thecorresponding transcription factor u
i
which resembles aconstant translation of mRNA to protein. The dynamics of the morphogen concentrationis then governed by a reaction-diffusion equation,
in which , , and d
j
, denote theconstant synthesis, degradation, and diffusion rate, respectively. In the following,we consider the anterior-posterior direction of the neural tube, which is large incomparison to the MHR and the expected diffusion length scale. Hence, no morphogenarrives at the border, i.e. we have zero morphogen concentration,v
j
(− L,t) = v
j
(L,t) = 0.The initial conditions are given by a two-dimensional vector of spatial functionsvj 0(x), j = 1,2,corresponding to the pattern at E8.5. Hence,
As the PDEs for v(t,x) are linearly coupled with the ODEsfor u(t,x), we can conclude that solutions exist for alltime points t ∈ [0,∞) [24]. Furthermore, the solution is unique, positive and in if the initial condition vector(u0,v0)T is positive andbounded [24]. This means that we obtain a solution which is continuous with respect totime and a square integrable function with respect to space and we can consider theasymptotic steady state of the system. Furthermore, the insight that the solutionsare square integrable in x ensures the convergence of finite-differencemethods, which are described in “Methods”. In the simulation, specialattention has to be paid to the non-differentiabilities occuring at least forOtx2 and Gbx2 in the stationary limit [12].
MHB model predicts stable steady state patterns
Biologically, one important characteristic of the GRN at the IsO is the activation ofthe genes in a precise spatio-temporal manner and the correct positioning of theirexpression domains [8]. Once the pattern is reached it has to be maintained even if the wholesystem is slightly perturbed for example by transcription noise. This implies thatthe model has to exhibit a stable stationary limit solution, which resembles thesteady state gene expression pattern that is observed at E10.5.
Due to the nonlinear coupling, the analytical calculation of the steady state isinfeasible for most components of the system. Only for Otx2 andGbx2 an analytical calculation of the steady state solution is possiblebecause they form a simple genetic toggle-switch system independent of the othercomponents. This system is well studied and we knew that it exhibits separatedexpression domains depending on the interaction strength of the two players [25, 26]. For simplicity and because we lacked detailed knowledge of theinteraction parameters of Otx2 and Gbx2, we assumed a symmetricparameter setup, i.e.
with α1 = α2 andβ1 = β2. Hereu1(t,x) denotes the expression ofOtx2 and u2(t,x) denotes theexpression of Gbx2. Based on the analytical calculation of the steady statewe found that the parameter values for which we obtained a steady state withseparated expression domains have to satisfy . Furthermore, the point x∗where both expression domains abut each other in the stationary limit is determinedby the initial conditions of Otx2, denoted byu01(x), and Gbx2, denoted byu02(x). Ifu01(x)>u02(x) forx<x∗ andu01(x)<u02(x) forx>x∗ then the boundary is placed atx∗.
In the following, the MHB is placed in the middle of the considered interval,x∗ = 0. Furthermore, similar to [10, 12], we set n = 2, k = 0.1, andα
i
= β
i
= 1for i = 1,…,8. The diffusion coefficients of Fgf8and Wnt1 were reduced compared to previous publications tod1 = d2 = 0.001.For the chosen diffusion parameters the values and are approximately zero, ensuring that the Dirichletboundary conditions have no significant influence on the dynamics or the steady stateof the systems, as assumed in the modeling process. For this setup we performed anextensive simulation study to determine the non-trivial steady states of the system(see “Methods” for details). If the initial transcription factor andmorphogen patters are unimodal, we find two stable, non-zero stationary limitsolutions. While the steady state 1 (Figure 2B, leftpanel) shows a maximal width for the expression domain of En andPax, for steady state 2 (Figure 2B, rightpanel) this width is minimal. The panels in Figure 2B showexemplary trajectories leading to the steady states. The simulation study indicates,that these are the only steady states reachable from initially unimodal expressiondomains.
For initially multimodal expression domains, we observed more complex expressionpattern, e.g., single spots of En or Pax expression. These arise asEn or Pax locally exceed the threshold of their mutualactivation, yielding isolated point regions with maximal En or Paxexpression, respectively. The spots are positioned in the region between the minimaland maximal width expression domain given by the steady states mentioned above. Wedid not consider the steady states with expression spots as they are not thought tobe of biological relevance, however this shows that the model is capable of producinga whole spectrum of steady states. Subsequently, we analyzed the stability of the twonon-trivial steady states to understand the temporal dynamics in their closeproximity. We found that both non-trivial steady states are stable. For both steadystates that we consider, the expression of Otx2, Gbx2 andWnt1 are the same, whereas for Fgf8, En andPax we observe a difference in the width of the individual expressiondomains.
Predicted steady states coincide qualitatively with the experimentalobservations
We compared our steady state solutions and their dependence on initial conditions andparameters to experimentally validated expression patterns. In our model the onlyfactor which determines the position of the MHB is the initial expression pattern ofOtx2 and Gbx2. This major role of Otx2 andGbx2 has already been shown in in vivo knock-out or knock-inexperiments, in which a change of Gbx2 and Otx2 expression domainsalso led to a change in the position of the MHB [27, 28]. Furthermore, both steady states show a restriction of the En andPax expression to the MHB region, which can also be observed invivo[8]. The steady state with the more restricted initial pattern(Figure 2B right panels and Figure 2C) shows a tighter expression domain. This tighter domain is due to theactivating interaction between En and Pax, whereas the other steadystate is dominated by the activating interaction of Fgf8 and Wnt1with En.
Concerning the expression domain of Fgf8 we found that in steady state 2 theexpression is restricted in posterior direction due to the regulation ofFgf8 by En and Pax. Only in this steady state a sharpboundary of the Fgf8 expression domain in posterior direction occurs(indicated by the dashed circles in Figure 2B,D). Thisagrees with the experimental observations that the expression of Fgf8 isrestricted to a ring at the caudal border of the MHB [29]. Furthermore, we concluded that the interactions which are not necessaryto maintain the expression pattern according to [10], are required to sharpen the expression domain of the morphogenFgf8. Those interactions are the activation of En byFgf8 and Wnt1 and the activation of Fgf8 by Enand Pax, which were not considered in [12].
Finally, we considered the expression pattern of Wnt1, which is the same inboth steady states. The IsO is a signaling center and one of its main tasks is togenerate a well-defined Wnt1 gradient. The gradient results from the sharplyrestricted and positioned expression domain of Wnt1. However, unfavorablesmooth interfaces of expression domains occur if expression is only regulated by amorphogen in our case Fgf8[30]. Depending on the diffusion coefficient of Fgf8 and theactivation of Wnt1 with Fgf8, the expression domain ofWnt1 becomes increasingly smooth. This disagreement with the experimentalresults where Wnt1 is expressed in a narrow ring at the rostral border ofthe MHB with a clearly visible boundary [8, 31–33]. For a detailed illustration and semi-quantitative assessment of theexpression domain of Wnt1 we refer to the EMAP eMouse Atlas Project [34] (http://www.emouseatlas.org/).
In the following section we will discuss possible post-transcriptional sharpeningmechanisms for the Wnt1 expression profile. As the more restrictedFgf8 expression pattern in steady state 2 is closer to the experimentalobservations [29], we will use this steady state in the following analysis, however, theresults are similar if steady state 1 is used. In particular the Wnt1expression patterns are alike for both steady states.
Sharpening of the Wnt1 expression by miRNA interactions
As our simulations showed no sharply delimited expression domain in the anteriordirection for the Wnt1 expression domain (see Figure 2B), there has to exist an unknown mechanism enforcing the sharpening ofthe Wnt1 profile over time, which is not considered in the model. In thiswork we consider post-transcriptional miRNA regulation as recent studies proved thatmiRNAs regulation is particular active at low mRNA copy numbers, which occur in ourmodel simulation at the Wnt1 expression domain boundary, and can induce geneexpression thresholds [30, 35]. This might results in a spatial sharpening of the target expressiondomain in the overlapping region [30, 35] via one of the following mechanisms:
-
i)
The miRNA binds transiently to the mRNA and only the mRNA is degraded (low degree of complementarity).
-
ii)
The mRNA-miRNA complex is degraded (high degree ofcomplementarity).
-
iii)
The mRNA-miRNA complex is degraded and unbound miRNAs are activelytransported between neighboring cells [36, 37].
The scenarios are illustrated in Figure 3A. It is alsoobserved that miRNAs do not lead to a degradation but to a translational inhibition,which would results in a complex accumulation and a low degradation rate in scenarioii) and iii). Mathematically, the model extensions are given by
with boundary conditions
and initial condition
Here, m(t,x) is the relative concentration of miRNA,c(t,x) is the relative concentration of mRNA-miRNAcomplex, and u4 is the relative Wnt1 mRNA level.β
m
is the degradation rate,d
m
is the diffusion coefficient of the miRNA andα
m
(x) is the space dependent productionprofile. Following the suggestions in [30] we used a sigmoidal shaped functionα
m
(x) = p1(tanh((l − x)/p2) + p4)to model the spatial dependence of the miRNA synthesis. The interaction strength ofmiRNA and Wnt1 mRNA is given by the binding rate κ, thedegradation rate of the resulting mRNA-miRNA complex is λ, and thedegree of miRNA recycling is denote by ξ.
Given this general model, the three scenario can be described by different parametersets. For scenario i), the turn-over parameters of the miRNA are set to zeroα
m
(x)≡0 andβ
m
= 0, assuming time-independentoverall miRNA level. Furthermore, miRNA is assumed to be completely recycled,ξ = 1, but not transported,d
m
= 0. In contrast, for scenario ii) andiii) ξ = 0,α
m
(x)≥0∀x ∈ U>0, andβ
m
>0. For these to parameters merely thediffusion coefficient is different, i.e. for ii)d
m
= 0 and for iii)d
m
>0. For all scenarios κ>0 andλ>0.
For this extended macroscopic model the existence, positivity and uniqueness is alsoguaranteed if and only if α
m
(x) is abounded, Lipschitz continuous function. This is the case as we assumedα
m
(x) to be a production profile withvalues in [0,1]. Hence, we could conclude that unique, positive solutions exist forthe extended model and we can simulate the model with the proposed algorithm. It isnot surprising that an overlap of the miRNA and Wnt1 profile leads to asharpening of the Wnt1 expression profile. However, we are especiallyinterested in the spatial shape an miRNA expression domain must have to sharpen theWnt1 boundary sufficiently as this can easily be compared to experimentalfindings and will lead to predictions for possible regulating miRNAs. In this contextwe defined a sharpening as reduction of Wnt1 in anterior direction and noreduction at the MHB, i.e. an increase in the second derivative with respect tox of the Wnt1 transcription level forx ∈ (0,L/2).
We simulated the three scenarios for different parameter sets(α
m
,d
m
,κ,λ)using the same initial conditions and parameters for the original model componentsand the profile functionα
m
(x) = 0.1(tanh((l − x)/0.3) + 1)with varying l (the profile is shown in Figure 3E). A representative simulation result is shown in Figure 3C). In addition, we studied the influence of the miRNA expressiondomain on the sharpening effect in the different scenarios. Therefore we varied theoverlap of the domains by increasing the profile function parameter l from 0(no-overlap) to 1 (constant production along the whole MHR). The results are shown inFigure 3D. We found that for scenario i) and ii) themiRNA level and the Wnt1 expression domain have to overlap in the regionwhere Wnt1 shows intermediate expression, i.e.l ∈ [0.4,0.5]. For both scenarios a complete overlap ofboth domains leads to a overall reduction of Wnt1 and we see no specificsharpening, especially for scenario i) a complete overlap led to a Wnt1knock-down state for the set of considered parameters. For scenario iii) the domainsmust only slightly overlap due to diffusion, if the domains strongly overlap themiRNA diffusing from the posterior direction leads to a blurring of the Wnt1domain at the MHB. We conclude that this scenario is not suitable for the MHB modelif we have a strongly overlapping miRNA domain. In the following, we focus onscenario i) for a low and ii) for a high degree of mRNA-miRNA complementarity [17–19].
miR-709 regulates Wnt1 mRNA expression in vitro
Up to this point, hypotheses about a possible regulation of Wnt1 by miRNAswas based on available knowledge about miRNA-mRNA interaction [30, 35] and simulation studies of the MHB model. To gain additional insight, asearch for experimental evidence of miRNAs possibly regulating the Wnt1 mRNAexpression was conducted. Therefore, we performed a miRNA target prediction andexperimentally validated the predicted miRNAs by determining their expression patternin the developing mouse embryo (with a special focus on the MHR) and their ability totarget the Wnt1 mRNA (3’UTR) in vitro.
We identified potential miRNAs targeting the Wnt1 mRNA using a combinationof several target prediction tools (see “Methods”). To reduce the numberof false-positive predictions, we post-processed the resulting list using logicfiltering. Therefore, we defined a logical state (ON/OFF) table for the three mouseembryonic stages, namely E8.5, E10.5 and E12.5, including Wnt1 and the otherfive IsO genes. E12.5 is considered as the IsO gene expression pattern, which hasrefined to its sharp domains at E10.5, is still maintained at E12.5. Assuming that afunctional miRNA targeting Wnt1 changes the amount of produced Wnt1 proteinand hence the expression state of genes downstream of Wnt1, we concluded that thosegenes change its expression state during miRNA expression. The target prediction wasfiltered for those miRNAs that target at minimum two “OFF genes” and no“ON genes” for each defined developmental stage. This resulted in a listof four miRNAs possibly regulating Wnt1 mRNA expression (Additional file1: Figure S1). From these possible candidate miRNAs(miR-705 and miR-709) were selected by ranking the interactionsaccording to the prediction scores provided by the distinct prediction tools.
To establish whether these two predicted miRNAs are expressed at the MHB in a patternthat is consistent with the model assumptions, their transcriptional profile in E10.5and E12.5 wild-type mouse embryos was determined. Therefore, we used a very sensitiveradioactive in situ hybridization method (for details see“Methods”) to detect the expression profile of these miRNAs along theentire anterior-posterior extent of the MHR on sagittal sections of theseembryos.
In the E10.5 mouse embryo, miR-705 is expressed only very weakly across theMHB and in the midbrain and rostralmost hindbrain, whereas miR-709 isexpressed strongly and uniformly across the MHB and in the midbrain and rostralhindbrain (see Additional file 2: Figure S2 and Additionalfile 2: Figure S3). Such a profile would lead to an overallreduction of Wnt1 mRNA in the model proposed but has no sharpening effect.In the E12.5 mouse embryo, by contrast, miR-705 was expressed strongly anduniformly in the entire MHR, including the midbrain, MHB and rostral hindbrain.However, we noted a graded expression of miR-709 across the MHB in theventral MHR (which is the region used to determine the expression profiles of theother IsO genes in the considered model) at this stage. Transcription ofmiR-709 at E12.5 was highest in the midbrain, declined towards the MHBand was lowest in the rostral hindbrain, the region of the hindbrain that is underthe influence of the IsO (see Figure 4A). This gradedmiR-709 expression profile became even more evident when we calculatedthe grayscale profile of miR-709 expression across the MHR, i.e. from theanterior end of the midbrain to the posterior end of the rostral hindbrain, as shownin Figure 4B (see “Methods”).
Following [30] we used the miRNA profile functionα
m
(x) = p1(tanh((l − x)/p2) + p3)and in order to analyze the sharpening effect of this profile we estimated theparametersp = (p1,p2,p3)and l (the resulting profile is shown in Figure 3E). We obtained a best-fit profile according to the grayscale profile(Figure 4B red profile). The profile is in accordancewith the modeling assumptions and we assumed that the expression pattern of the IsOgenes is already stably established at E10.5 and this particular miRNA profile is notestablished before E12.5. This is evidence that the miRNA regulation of Wnt1in the model is not a regulatory interaction necessary to establish the pattern, butrather acts as a fine tuning mechanism to reduce the range of cells which have anintermediate Wnt1 expression. To verify this effect the model was simulatedusing the E10.5 expression pattern as initial condition and the estimated profile formiRNA production. We used scenario ii) for the simulation (see Figure 4C and D), because scenario i) mostly affects the repression ofWnt1 translation and we only have evidence for degradation with theperformed experiments. We also neglect iii) as we had no experimental evidence for adiffusion of miR-709 in the neural tube. In the simulations, a clearlyvisible sharpening effect could be observed, especially if we increase theFgf8 diffusion coefficient (see Figure 4C andD).
In contrast to miR-709, did not show a refined and graded expression aroundthe ventral MHB in the E12.5 mouse embryo (see Additional file 3: Figure S3). These results indicated that only miR-709 isexpressed within the MHR and around the MHB in a pattern as predicted by the modeland consistent with a possible regulatory function of miR-709 onWnt1 expression at the MHB.
To determine whether miR-709 and miR-705 can indeed regulate theexpression of Wnt1 in a cellular context, we used the so-called luciferase“sensor assays”. In this experimental setup, a fragment of the Wnt13’UTR containing the predicted miR-709 and miR-705binding sites (BS) was cloned at the 3’ end of a sequence encoding the Fireflyluciferase protein. This constitutes the so-called “sensor vector”. Theluciferase protein transfected with the sensor vector are indirectly measured by abioluminescence reaction resulting in the emission of light, and the intensity of theemitted light is directly proportional to the luciferase protein concentration in thecells. Co-transduction of the cells with the Wnt1 3’UTR sensor vectorand miR-709 or miR-705 should result in a reduction of luciferaseprotein levels, relative to an “empty” control vector (that does notcontain any Wnt1 3’UTR sequences), if and only if these miRNAs bind totheir target sites within the Wnt1 3’UTR and post-transcriptionallyinhibit the expression of luciferase protein in these cells. Indeed, co-transductionof HEK-293T cells with the Wnt1 3’UTR sensor vector andmiR-709 resulted in an approx. 23% reduction of luciferasebioluminiscence in these cells (see Additional file 3:Figure S2), whereas co-transduction of HEK-293T cells with the Wnt13’UTR sensor vector and miR-705 had no significant effect (seeAdditional file 3: Figure S3). This result indicated thatmiR-709, but not miR-705, indeed targets theWnt13’UTR. Next, we determined whether the post-transcriptionalregulation of Wnt1 expression by miR-709 is indeed mediated by thepredicted miR-709 BS within the Wnt1 3’UTR. Therefore, thepredicted BS sequences were changed such that they could not be bound anymore bymiR-709 (see Additional file 4: Table S1) andconsequently the expression levels of luciferase protein should not be affectedrelative to the control (“empty”) sensor vector. Mutagenesis of the twopredicted miR-709 BS within the Wnt1 3’UTR in fact abolishedthe negative regulation of luciferase expression from the sensor vector bymiR-709. This result indicated that the negative post-transcriptionalregulation of the Wnt1 mRNA by miR-709 is in fact mediated by thetwo predicted miR-709 BS within the Wnt1 3’UTR.