A robust mathematical model of adaxial–abaxial patterning

Abstract

Biological development results from intricate and dynamic interactions between members of gene regulatory networks. This is exemplified by the production of flat leaf architecture. Leaves flatten by driving growth along the boundary between their adaxial (top) and abaxial (bottom) domains. These domains are generated by interactions between a complex network of transcription factors and small RNAs. Despite its complexity, flat leaf production is robust to genetic and environmental noise. To identify factors contributing to this robustness, we mathematically modelled the determinants and interactions that pattern the adaxial–abaxial axis in leaves of Arabidopsis thaliana. Model parameters were estimated almost exclusively using experimental data. Our model recapitulates observations of adaxial–abaxial patterning and small RNA-target interactions. Positioning of the adaxial–abaxial boundary is stable across a wide range of small RNA source values and is highly robust to noise in the model. The successful application of our one-dimensional spatial model will enable higher-dimension modelling of the complex and mechanistically challenging process of flat leaf production.

1. INTRODUCTION

A wide variety of complex and challenging processes unfold in a highly reproducible (or robust) manner. In biology, this phenomenon is exemplified by developmental patterning, which reliably generates complex structures, even in the face of genetic and environmental noise. For instance, non-linear relationships between inputs and outputs are central to robust systems (Félix and Barkoulas 2015; Green et al. 2017), and some mechanisms facilitating this have been identified (Rutherford and Lindquist 1998; Félix and Wagner 2008; Albergante et al. 2014). However, molecular mechanisms underlying developmental reproducibility remain largely unknown.

Plant development is a particularly attractive system to tackle these questions. As plants are fixed to one location, they must be highly sensitive to environmental stimuli yet still robustly translate information into discrete, context-dependent outcomes. This is illustrated by the production of flat leaf architecture. Plants carry a stem cell niche from which leaves and flowers continuously emerge throughout their life cycle (Maughan et al. 2006). Further, leaves do not start out flat but instead emerge from the stem cell niche as radially symmetric bulges. These bulges then flatten by driving growth along two of its axes to create final structures that are long and wide, but only a few cell layers thick (Husbands et al. 2015). Despite the difficulty of this process, plants consistently form leaves with nearly identical geometries. Flat leaf architecture is therefore an excellent system to provide insights into mechanisms that lend robustness to biological processes.

Leaves differentiate along multiple developmental axes, including top to bottom (adaxial to abaxial). The leaves of the model plant Arabidopsis thaliana are typically comprised of six cell layers (two adaxial and four abaxial), and each domain contains distinct cell types with specialized roles in development (Fig. 1A and B) (Waites and Hudson 1995; Eshed et al. 2001; Engstrom et al. 2004). For example, adaxial cells are rectangular and densely packed to maximize light capture, while abaxial cells are irregularly shaped and loosely packed to facilitate gas exchange. This spatial arrangement optimizes photosynthesis while minimizing water loss (Engstrom et al. 2004). The adaxial–abaxial axis is also critical to flat leaf production, as leaves use the boundary between their adaxial and abaxial sides to orient directional growth. Even minor defects in the coordination of this growth compromise leaf flatness, leading to severe fitness consequences (Lang et al. 2004; Zhang et al. 2009; Nakata et al. 2012). Robust flat leaf production thus depends on a properly patterned adaxial–abaxial axis, which in turn is a reflection of robust outputs from the underlying network.

Figure 1.

Adaxial and abaxial layers are molecularly and morphologically distinct. (A, B) Adaxial (AD) and abaxial (AB) domains are morphologically distinct. Arabidopsis leaves usually possess two layers of adaxial cells and four layers of abaxial cells. (C) The adaxial–abaxial axis is patterned by a complex and redundant network of transcription factors and small RNAs (see main text). T-bars denote repressive relationships. Direct interactions (solid lines) and indirect interactions (dotted lines) are noted. Barbell between ARF and KAN denotes protein–protein interaction. (D) Cell layers (C₁ through C₆) of the adaxial–abaxial axis constituting the fixed domain of the model. Shown are sites of accumulation for mobile small RNAs and cell-autonomous transcription factors. Triangles extending from miR166 and tasiARF transcription domains illustrate small RNA diffusion gradients.

Open in new tab Download slide

Mathematical modelling is a powerful approach to study complex developmental systems (Tomlin and Axelrod 2007; Sharpe 2017). In plants, mathematical modelling has been used to examined aspects of development ranging from hormone signalling to cell-fate specification to the acquisition of complex floral organ morphologies (Tameshige et al. 2013; Muraro et al. 2014; Miyashima et al. 2019). Recent work has also modelled the role of the adaxial–abaxial axis in the production of flat leaves and the invaginated feeding structures of carnivorous plants (Fukushima et al. 2015; Hayakawa et al. 2016; Whitewoods et al. 2020). These models highlight the importance of the adaxial–abaxial axis in orienting growth and division, and suggest the division plane as a key determinant in producing flattened outgrowth. However, they do not include interactions between the extensive network of adaxial–abaxial polarity determinants, and it remains unclear how the division plane itself is established and maintained. Further, mechanisms that might explain the robustness of these complex biological processes are not discussed.

To address this knowledge gap, we developed a one-dimensional spatial model of adaxial–abaxial patterning. Importantly, parameters were estimated almost exclusively using experimental data, and the resulting model is highly robust to noise. We find that the interaction between CLASS III HOMEODOMAIN-LEUCINE ZIPPER (HD-ZIPIII) transcription factors and miR166 is critical for robust specification of the adaxial–abaxial boundary. Taken together, our model suggests possible mechanisms underlying the robustness of flat leaf production and will inform future modelling of the numerous complex structures generated by adaxial–abaxial patterning.

2. METHODS

2.1 Identification and interaction of important substances

The adaxial–abaxial axis is patterned by a complex network of transcription factors and small RNAs (Husbands et al. 2009) (Fig. 1C and D). This network is characterized by extensive redundancy, as cell fate is often specified by multiple members of a given gene family, as well as by other distinct gene families. Dominant contributors to adaxial fate include the HD-ZIPIII family of deeply conserved homeodomain-containing transcription factors (McConnell et al. 2001; Juarez et al. 2004; Romani et al. 2018). Adaxial fate is also specified by the LOB DOMAIN transcription factor ASYMMETRIC LEAVES2 (AS2), which forms a complex with the Myb transcription factor ASYMMETRIC LEAVES1 (AS1) (Lin et al. 2003; Husbands et al. 2007; Iwakawa et al. 2007). Finally, the tasiARF trans-acting siRNAs produced from the TAS3A locus also play an important role in specifying adaxial fate (Allen et al. 2005; Nogueira et al. 2007). Abaxial fate is similarly specified by multiple transcription factors and small RNAs. Transcription factors include members of the KANADI (KAN) family, as well as AUXIN RESPONSE FACTOR3 and 4 (ARF3/ARF4) (Eshed et al. 2001; Kerstetter et al. 2001; Pekker et al. 2005), and physical interaction between KAN and ARF transcription factors is critical for their function (Kelley et al. 2012). The miR165/miR166 family of miRNAs also specifies abaxial fate (McConnell et al. 2001; Juarez et al. 2004; Mallory et al. 2004; Alvarez et al. 2006). Lastly, the WUSCHEL-RELATED HOMEOBOX (WOX) transcription factor is not a canonical polarity determinant per se, but rather affects the activity of both adaxial and abaxial determinants (Nakata et al. 2012).

The adaxial–abaxial network resolves into a system in which adaxial and abaxial cell fates are mutually exclusive (Fig. 1D). This is accomplished by extensive mutual antagonism between adaxial and abaxial polarity determinants. For instance, HD-ZIPIII proteins are transcribed throughout the leaf primordia but only accumulate in the adaxial domain (McConnell et al. 2001; Skopelitis et al. 2017). This patterning is derived from the activity of miR166, which is transcribed in the abaxial epidermis and forms a mobility-derived gradient that eliminates HD-ZIPIII transcripts from the abaxial side of the leaf (Juarez et al. 2004; Kidner and Martienssen 2004; Nogueira et al. 2007; Zhu et al. 2011). HD-ZIPIII proteins also indirectly repress miR166 activity by upregulating AGO10. This catalytically dead Argonaute effector protein sequesters miR166 into non-functional complexes (Zhu et al. 2011; Brandt et al. 2013). Loss of miR166-mediated cleavage of HD-ZIPIII transcripts results in severely adaxialized phenotypes (McConnell et al. 2001; Mallory et al. 2004; Alvarez et al. 2006), demonstrating the importance of the HD-ZIPIII–MIR166 relationship to adaxial–abaxial patterning.

ARF3 is similarly patterned by tasiARFs, as uniformly transcribed ARF3 is restricted to the four abaxial cell layers by a gradient of tasiARFs produced from the adaxial side of leaf primordia (Chitwood et al. 2009). KAN genes are transcribed in the bottom two abaxial layers and function by directly repressing AS2 and indirectly repressing WOX (Wu et al. 2008; Nakata et al. 2012). Repression of AS2 is also redundantly mediated by WOX transcription factors, which directly restrict transcription of AS2 to the adaxial epidermis (Nakata et al. 2012; Zhang et al. 2014). In complex with AS1, AS2 directly represses multiple abaxial determinants, including miR166 and ARF, using a variety of molecular mechanisms (Husbands et al. 2015). These substances and interactions comprise the foundation of our mathematical model.

2.2 Equations used in model

In constructing our one-dimensional spatial model of Arabidopsis leaf development, we describe quantities of HD-ZIPIII, AS2, tasiARF, miR166, ARF, KAN and WOX gene classes in a single column of six leaf cells (Fig. 1B). Quantities are measured in mRNA transcripts per million (TPM), which reports normalized transcript abundance in each cell after controlling for sequencing depth. We denote these values using the functions H, S, T, M, R, K and W, respectively. Note that both transcription factor and small RNA gene classes require additional steps before they can contribute to adaxial–abaxial patterning. Specifically, mRNAs of transcription factors are translated into proteins, while mature small RNAs are produced through a complex multistep process (Yu et al. 2017). Incorporating these additional steps would significantly complicate the model and require us to make a number of unsupported assumptions. As such, the value of each function is taken as a direct proxy of activity.

We first selected a biologically realistic spatial domain on which to model these interactions. The length of the adaxial–abaxial axis in Arabidopsis can vary dramatically depending on the developmental stage (Nakata et al. 2012; Caggiano et al. 2017). We assumed a length of 20 μm from the adaxial epidermis to the abaxial epidermis, approximating the length of an early leaf primordium (Tameshige et al. 2013). Therefore, the spatial domain is taken to be [0, 20]. The intervals C_i = [20(i − 1)/6, 20i/6] for $i = 1, \dots, 6$ represent the six leaf cells. C₁ represents the adaxial epidermis, and C₆ represents the abaxial epidermis. We refer to the corresponding cells as cell C₁ through cell C₆.

The quantities of gene classes in each of the six cells are governed by the following ordinary differential equation (ODE) system:

\begin{array}{l} \frac{\partial H_{j}}{\partial t} = σ_{H} - d_{H} H_{j} - k_{H M} h_{H M}^{j} H_{j} M_{j} \end{array}

(1)

\begin{array}{l} \frac{\partial S_{j}}{\partial t} = \frac{σ_{S}}{1 + k_{W S} W + k_{K R S} K_{j} R_{j}} {x \in C_{1} \cup C_{6}} - d_{S} S_{j} \end{array}

(2)

\begin{array}{l} \frac{\partial T_{j}}{\partial t} = D_{T} Δ_{e} T_{j} + σ_{T} {x \in A d} - d_{T} T_{j} - k_{T R} h_{T R}^{j} T_{j} R_{j} \end{array}

(3)

\begin{array}{l} \frac{\partial M_{j}}{\partial t} = D_{M} Δ_{e} M_{j} + \frac{σ_{M}}{1 + k_{S M} S_{j}} {x \in C_{6}} - d_{M} M_{j} - k_{H M} h_{H M}^{j} H_{j} M_{j} \end{array}

(4)

\begin{array}{l} \frac{\partial R_{j}}{\partial t} = \frac{σ_{R}}{1 + k_{S R} S_{j}} - d_{R} R_{j} - k_{T R} h_{T R}^{j} T_{j} R_{j} \end{array}

(5)

\begin{array}{l} \frac{\partial K_{j}}{\partial t} = σ_{K} {x \in A b} {x \in C_{5} \cup C_{6}} - d_{K} K_{j} \end{array}

(6)

\begin{array}{l} \frac{\partial W_{j}}{\partial t} = \frac{σ_{W}}{1 + k_{K W} K_{j}} {x \in A b} {t > t_{1}} - d_{W} W_{j} \end{array}

(7)

In the above equations, P_j refers to the TPM of substance P in cell C_j for $j = 1, \dots, 6$ ⁠. Given a Boolean expression A, the function {A} returns the value 1 when A is true and returns 0 otherwise. The abaxial region is defined to be Ab = {M > H}, and the adaxial region is defined to be Ad = 1 − Ab. Activation time for W is given by t₁ > 0.

The system also uses two Hill functions which are given as follows:

\begin{array}{l} h_{H M}^{j} = \frac{M_{j}^{n_{M}}}{H_{j}^{n_{M}} + M_{j}^{n_{M}}}, h_{T R}^{j} = \frac{T_{j}^{n_{T}}}{R_{j}^{n_{T}} + T_{j}^{n_{T}}} . \end{array}

Finally, Δ _e is a Laplacian operator which describes diffusion into adjacent cells (irrespective of mechanism) with Neumann boundary conditions. For example, $Δ_{e} T_{1} = (T_{2} - T_{1}) / Δ x^{2}$ ⁠, $Δ_{e} T_{j} = (- T_{j - 1} + 2 T_{j} - T_{j + 1}) / Δ x^{2}$ for $2 \leq j \leq 5$ and $Δ_{e} T_{6} = (T_{5} - T_{6}) / Δ x^{2}$ ⁠. Here, Δx = 10/3 μm, which is the length of each cell.

2.3 Terms used in equations

Our model of adaxial–abaxial patterning takes into account the following known or assumed properties of the biological system. First, adaxial and abaxial determinants are expressed in specific locations at specific times. Second, transcription factors in the network often function by repressing transcription of opposing determinants. Third, interaction between small RNAs and their targets is repressive and destructive, leading to the turnover of both molecules. Fourth, small RNAs are mobile and form gradients across the leaf with the highest concentration located at their site of biogenesis. Fifth, these small RNA gradients generate binary readouts of their targets in the leaf, reminiscent of protein morphogenic behavior (Chitwood et al. 2009; Skopelitis et al. 2017). Finally, proteins and small RNAs are naturally turned over independent of interactions with other determinants. These properties are modelled primarily using standard modelling techniques for enzyme kinetics in biological systems (Murray 2007; Chou and Friedman 2016).

Expression of adaxial–abaxial determinants is governed by source terms which contain a σ parameter. The value of this parameter corresponds to the rate of transcription in the absence of repression by other substances, with a value of zero assigned for cell layers not expressing the determinant. From established spatial expression patterns, we set H and R to be transcribed in all cell layers (McConnell et al. 2001; Chitwood et al. 2009; Skopelitis et al. 2017), S transcribed in layers 1 and 6 (epidermis) (Nakata et al. 2012; Zhang et al. 2014; Husbands et al. 2015), T transcribed in both adaxial cell layers (Chitwood et al. 2009), M transcribed in layer 6 (abaxial epidermis) (Husbands et al. 2015), K transcribed in layers 5 and 6 (two abaxial-most layers) (Wu et al. 2008; Nakata et al. 2012) and W transcribed in all four abaxial cell layers (Nakata et al. 2012). In addition, genetic and reporter assays suggest HD-ZIPIII, ARF, miR166 and KAN are present at the earliest stage of leaf development, termed the P0 (Heisler et al. 2005; Chitwood et al. 2009). WOX on the other hand does not accumulate until the P3 stage of leaf development, which begins ~72 h post-initiation (Nakata et al. 2012). We therefore include a factor of {t > t₁} in the source term for W, where t₁ = 72 h is the activation time for transcription of W. Specific transcription activation times for the rest of the genes are not known. However, our model reaches the same steady state regardless of substance activation times. As such, we assumed that all substances except for W are transcribed at the start of the model.

Interactions between adaxial–abaxial determinants are a key feature of the network (Husbands et al. 2015). Many transcription factors in the network affect the steady-state levels of their targets by direct repression of transcription. As such, repression of proteins by other proteins is incorporated in the denominators of source terms. Note that we include the term $k_{K R S} K R$ in the denominator of the S source term because KAN proteins physically interact with ARF proteins (Kelley et al. 2012). By contrast, small RNAs and their target transcripts physically interact and are turned over. These interactions are thus described by additional functions which use k parameters to denote repressive relationships. These terms are included in both equations as the interaction is assumed to destroy both substances at the same time (Levine et al. 2007; Muraro et al. 2014). For this reason, the equations for H_j (Equation (1)) and M_j (Equation (4)) both include the term $k_{H M} h_{H M}^{j} H_{j} M_{j}$ ⁠. Similarly, the equations for both T_j (Equation (3)) and R_j (Equation (5)) include the term $k_{T R} h_{T R}^{j} T_{j} R_{j}$ ⁠.

Transcription factors in the adaxial–abaxial network, represented by H, S, R, K and W, accumulate in the nucleus, and currently no experimental evidence suggests intracellular movement of small RNAs contributes to adaxial–abaxial patterning. As such, each cell is represented using a single point. Further, as transcription factors in the network function cell autonomously, they were not assigned extracellular diffusion operators. By contrast, the two small RNAs in the network, represented by T and M, are mobile (Chitwood et al. 2009; Skopelitis et al. 2017). As such, they were modelled using Δ _e diffusion operators to reflect this mobility.

Strikingly, small RNA gradients in the leaf do not lead to inverse gradients of their targets. Instead, these graded inputs produce binary on–off readouts (Chitwood et al. 2009; Skopelitis et al. 2017). This implies that in cells containing low levels of small RNAs and high levels of their targets, these substances can coexist without interacting. By contrast, in cells containing high levels of small RNAs, interaction is triggered between these substances leading to mutual repression. This sigmoidal or switch-like behavior is reminiscent of morphogens in animal systems, which are typically modelled using Hill functions (Weiss 1997; Park et al. 2019). For this reason, we included Hill functions in the interaction terms between small RNAs and their targets.

Finally, proteins and small RNAs are assumed to naturally turn over at particular rates which are modelled using d parameters. Note this is considered an inherent property of each substance and not a result of interactions with other substances in the model.

2.4 Parameters

In all model simulations, unless stated otherwise, we use the parameter values stated in Table 1. Units are given in the same table. The model is implemented using a forward Euler method in time [see Supporting Information—File S1], which can be implemented using standard techniques (LeVeque 2007). Our simulations used a time step of Δt = 0.01 h, which is small enough to ensure stability of our numerical method.

Table 1.

Open in new tab

Parameter values

Diffusion parameters	Value	Studied range	Reference
$D_{T}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
$D_{M}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
Degradation parameters	Value	Studied range	Reference
$d_{H}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{S}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{T}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{M}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{R}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{K}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{W}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
Source parameters	Value	Studied range	Reference
$σ_{H}$	25.8 h⁻¹		Cortijo et al. (2019)
$σ_{S}$	3.6 h⁻¹		Cortijo et al. (2019)
$σ_{T}$	8.635 h⁻¹	6.24–11.03 h⁻¹
$σ_{M}$	223.2 h⁻¹	189.3–257.1 h⁻¹
$σ_{R}$	6 h⁻¹		Cortijo et al. (2019)
$σ_{K}$	8.4 h⁻¹		Cortijo et al. (2019)
$σ_{W}$	3 h⁻¹		Cortijo et al. (2019)
Interaction parameters	Value	Studied range	Reference
$k_{H M}$	1 TPM⁻¹ h⁻¹
$k_{T R}$	1 TPM⁻¹ h⁻¹
$k_{K W}$	2 TPM⁻¹ h⁻¹
$k_{S M}$	2 TPM⁻¹ h⁻¹
$k_{S R}$	1 TPM⁻¹ h⁻¹
$k_{W S}$	2 TPM⁻¹ h⁻¹
$k_{K R S}$	2 TPM⁻¹ h⁻¹
Hill coefficients	Value	Studied range	Reference
$n_{M}$	6
$n_{T}$	6
Activation time	Value	Studied range	Reference
$t_{1}$	72 h		Nakata et al. (2012)

Diffusion parameters	Value	Studied range	Reference
$D_{T}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
$D_{M}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
Degradation parameters	Value	Studied range	Reference
$d_{H}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{S}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{T}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{M}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{R}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{K}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{W}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
Source parameters	Value	Studied range	Reference
$σ_{H}$	25.8 h⁻¹		Cortijo et al. (2019)
$σ_{S}$	3.6 h⁻¹		Cortijo et al. (2019)
$σ_{T}$	8.635 h⁻¹	6.24–11.03 h⁻¹
$σ_{M}$	223.2 h⁻¹	189.3–257.1 h⁻¹
$σ_{R}$	6 h⁻¹		Cortijo et al. (2019)
$σ_{K}$	8.4 h⁻¹		Cortijo et al. (2019)
$σ_{W}$	3 h⁻¹		Cortijo et al. (2019)
Interaction parameters	Value	Studied range	Reference
$k_{H M}$	1 TPM⁻¹ h⁻¹
$k_{T R}$	1 TPM⁻¹ h⁻¹
$k_{K W}$	2 TPM⁻¹ h⁻¹
$k_{S M}$	2 TPM⁻¹ h⁻¹
$k_{S R}$	1 TPM⁻¹ h⁻¹
$k_{W S}$	2 TPM⁻¹ h⁻¹
$k_{K R S}$	2 TPM⁻¹ h⁻¹
Hill coefficients	Value	Studied range	Reference
$n_{M}$	6
$n_{T}$	6
Activation time	Value	Studied range	Reference
$t_{1}$	72 h		Nakata et al. (2012)

Values for all model parameters. Measured parameter ranges and corresponding references are provided when available.

Table 1.

Open in new tab

Parameter values

Diffusion parameters	Value	Studied range	Reference
$D_{T}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
$D_{M}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
Degradation parameters	Value	Studied range	Reference
$d_{H}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{S}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{T}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{M}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{R}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{K}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{W}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
Source parameters	Value	Studied range	Reference
$σ_{H}$	25.8 h⁻¹		Cortijo et al. (2019)
$σ_{S}$	3.6 h⁻¹		Cortijo et al. (2019)
$σ_{T}$	8.635 h⁻¹	6.24–11.03 h⁻¹
$σ_{M}$	223.2 h⁻¹	189.3–257.1 h⁻¹
$σ_{R}$	6 h⁻¹		Cortijo et al. (2019)
$σ_{K}$	8.4 h⁻¹		Cortijo et al. (2019)
$σ_{W}$	3 h⁻¹		Cortijo et al. (2019)
Interaction parameters	Value	Studied range	Reference
$k_{H M}$	1 TPM⁻¹ h⁻¹
$k_{T R}$	1 TPM⁻¹ h⁻¹
$k_{K W}$	2 TPM⁻¹ h⁻¹
$k_{S M}$	2 TPM⁻¹ h⁻¹
$k_{S R}$	1 TPM⁻¹ h⁻¹
$k_{W S}$	2 TPM⁻¹ h⁻¹
$k_{K R S}$	2 TPM⁻¹ h⁻¹
Hill coefficients	Value	Studied range	Reference
$n_{M}$	6
$n_{T}$	6
Activation time	Value	Studied range	Reference
$t_{1}$	72 h		Nakata et al. (2012)

Diffusion parameters	Value	Studied range	Reference
$D_{T}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
$D_{M}$	50 μm² h⁻¹	9–360 μm² h⁻¹	Valdez-Taubas and Pelham (2003)
Degradation parameters	Value	Studied range	Reference
$d_{H}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{S}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{T}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{M}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{R}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{K}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
$d_{W}$	0.6 h⁻¹	0–1.3 h⁻¹	Sorenson et al. (2018)
Source parameters	Value	Studied range	Reference
$σ_{H}$	25.8 h⁻¹		Cortijo et al. (2019)
$σ_{S}$	3.6 h⁻¹		Cortijo et al. (2019)
$σ_{T}$	8.635 h⁻¹	6.24–11.03 h⁻¹
$σ_{M}$	223.2 h⁻¹	189.3–257.1 h⁻¹
$σ_{R}$	6 h⁻¹		Cortijo et al. (2019)
$σ_{K}$	8.4 h⁻¹		Cortijo et al. (2019)
$σ_{W}$	3 h⁻¹		Cortijo et al. (2019)
Interaction parameters	Value	Studied range	Reference
$k_{H M}$	1 TPM⁻¹ h⁻¹
$k_{T R}$	1 TPM⁻¹ h⁻¹
$k_{K W}$	2 TPM⁻¹ h⁻¹
$k_{S M}$	2 TPM⁻¹ h⁻¹
$k_{S R}$	1 TPM⁻¹ h⁻¹
$k_{W S}$	2 TPM⁻¹ h⁻¹
$k_{K R S}$	2 TPM⁻¹ h⁻¹
Hill coefficients	Value	Studied range	Reference
$n_{M}$	6
$n_{T}$	6
Activation time	Value	Studied range	Reference
$t_{1}$	72 h		Nakata et al. (2012)

Values for all model parameters. Measured parameter ranges and corresponding references are provided when available.

Wherever possible, we used primary literature to identify experimentally derived or biologically relevant values to inform estimates of our parameters. This is not possible for small RNA diffusion parameters, as these values have not been experimentally determined. However, we have previously constructed models faithfully replicating biological phenomena in yeast (Chou et al. 2008, 2011; Renardy et al. 2018). Importantly, these models accurately simulate biology using experimentally derived diffusion parameters, which lie between 9 and 360 μm² h⁻¹ (0.0025 and 0.1 μm² s⁻¹) (Valdez-Taubas and Pelham 2003). We thus tested multiple parameters in this range, and found that diffusion parameters of 50 μm² h⁻¹ were sufficient to generate a biologically realistic model.

To estimate degradation parameters we used experimentally determined degradation rates for 18 000 Arabidopsis mRNAs (Sorenson et al. 2018). This data set has an average of ~0.6 h⁻¹ and standard deviation of ~0.7 h⁻¹. As such, we set all degradation parameters to 0.6 h⁻¹ with a studied range of 0–1.3 h⁻¹.

To compute source parameters, suppose the quantity of a substance G is modelled using the ODE

\begin{array}{l} G^{'} (t) = σ_{G} - d_{G} G (t) \end{array}

(8)

where $σ_{G}, d_{G} > 0$ are constants. Given initial condition $G (0) = G_{0}$ ⁠, the exact solution is given by

\begin{array}{l} G (t) = (G_{0} - σ_{G} / d_{G}) e^{- d_{G} t} + σ_{G} / d_{G} . \end{array}

(9)

This equation has a global steady state of $σ_{G} / d_{G}$ ⁠. Therefore, using estimated steady-state values and degradation parameters, we can compute the σ parameters. To estimate protein steady states, we use data from a study measuring global transcript abundance of individual Arabidopsis seedlings at various time points (Cortijo et al. 2019). This provides the relative abundance in all cell types of each determinant (measured in TPM). To compute average TPM per cell, we normalized for the size of the expression domain and used the resulting values rounded to nearest integers as our steady states [see Supporting Information—File S2]. Source parameters were then obtained by multiplying these steady-state values by their degradation parameters. For example, the total TPM value for the HD-ZIPIII genes is ~86. Since HD-ZIPIII is restricted to the two cells of the adaxial domain, the TPM of HD-ZIPIII per cell in its accumulation domain is 43, and $σ_{H} = 43 d_{H} = 25.8$ h⁻¹.

As the MIR166 genes that code for miR166 were not represented in the data set (Cortijo et al. 2019), we determined small RNA source values using an approach based on outputs of the model itself. To determine $σ_{M}$ ⁠, we fixed the other model parameters and mapped the range of $σ_{M}$ values which resulted in a biologically realistic output of adaxial–abaxial patterning. To be considered biologically realistic, this must satisfy two conditions. First, the adaxial–abaxial boundary must generate precisely two adaxial cell layers. Second, there must be near-equivalent levels of HD-ZIPIII transcripts in these two adaxial layers (i.e. no >10 % deviation), in accordance with experimental evidence (Skopelitis et al. 2017) (discussed in Results). The maximal interval J₂ such that the endpoints of J₂ are accurate to one decimal place and the model satisfies the above two conditions for all $σ_{M} \in J_{2}$ is given by J₂ = [189.3, 257.1]. For our specific parameter value, we chose $σ_{M} = 223.2$ h⁻¹, which is the midpoint of J₂. A similar strategy was employed for T, which found an interval of [6.24, 11.03]. The midpoint of this interval, $σ_{T} = 8.635$ h⁻¹, was then selected as our parameter value.

$n_{M}$ and $n_{T}$ are Hill coefficients that describe the behavior of miR166 and tasiARFs, respectively. As experimental data are not available for either coefficient, we again used model outputs to inform our selections. We tested multiple values for each parameter and scored the frequency at which the resulting models displayed biologically realistic behavior (defined in previous section). The adaxial–abaxial boundary is thought to be established by the activities of miR166 and HD-ZIPIII (Husbands et al. 2009; Skopelitis et al. 2017). As such, we used H and M to evaluate possible $n_{M}$ parameter values (Fig. 2). For each $n_{M}$ value, we used the smallest integer value of $σ_{M}$ such that the model results in two adaxial cells. In the absence of the Hill function (⁠ $n_{M} = n_{T} = 0$ ⁠), cell C₂ contains notably smaller values of H than cell C₁, and cell C₃ is of uncertain fate (Fig. 2A). Neither observation matches experimental data. At $n_{M} = 2$ ⁠, cells C₂ and C₁ have near-equivalent values of H; however, cell C₃ still does not definitively acquire either adaxial or abaxial fate (Fig. 2B). Only at $n_{M} = 4$ do we observe equivalent values of H in cells C₂ and C₁, and cell C₃ can be clearly distinguished as possessing abaxial fate (Fig. 2C). By $n_{M} = 6$ ⁠, values of H in cell C₃ drop to nearly undetectable levels (Fig. 2D). $n_{M} = 6$ was therefore selected as our Hill coefficient value.

Figure 2.

Hill coefficients are required to recapitulate observed HD-ZIPIII–miR166 accumulation patterns. (A) $n_{M} = 0$ yields uneven levels of HD-ZIPIII (red) in adaxial cells and nearly identical levels of HD-ZIPIII and miR166 (green) in cell C₃. As $n_{M}$ increases to (B) 2, (C) 4 and (D) 6, HD-ZIPIII levels in adaxial cells equilibrate, and cell C₃ accumulates higher levels of miR166 than HD-ZIPIII.

Open in new tab Download slide

Although we modelled different molecular interactions, our findings are consistent with known properties of Hill functions. For example, Hill coefficients >8 are considered biologically unrealistic (Gonze and Abou-Jaoudé 2013). Our model recapitulates observed HD-ZIPIII–miR166 behavior using Hill coefficients well below that value. Perhaps more importantly, our independently obtained Hill coefficient closely matches experimentally measured values in other systems. For example, Hill functions used to model morphogen interactions in Drosophila melanogaster found Hill coefficients of 6.2 ± 1.1 and 5.2 ± 0.3 for hunchback and WTHbP2, respectively (Park et al. 2019). When applied to our adaxial–abaxial patterning model, both of these Hill coefficient ranges result in realistic behavior (Fig. 2D), further supporting the idea that the Hill coefficients necessary for our model are biologically realistic.

Wherever possible, parameters were selected using available quantitative evidence. However, this was not possible for some parameters. The values of these parameters were instead chosen such that model output qualitatively matched experimental evidence. Thus, while some model parameters were selected qualitatively and others quantitatively, in all cases, the choice of parameter values was directly informed by experimental evidence.

3. RESULTS

3.1 Adaxial–abaxial patterning can be mathematically modelled

Our model accurately recapitulates observations of adaxial–abaxial patterning in wild-type leaves (Fig. 3; see Supporting Information—File S3) (Husbands et al. 2009, 2015). For instance, at the steady state of t = 96 h, the activity of the uniformly produced H becomes localized to two adaxial cell layers by M. R, which is also produced throughout the primordium, is localized to four abaxial layers via T activity. Thus, the final shape functions of H and R resemble those of HD-ZIPIII and ARF3, respectively. The accumulation patterns of M and T also resemble those of their wild-type counterparts. miR166 is expressed at high levels, is found at much higher concentration in abaxial cells than adaxial cells and sharply drops off between cell layers (Skopelitis et al. 2017). By contrast, tasiARF levels are low, concentrated in adaxial cells and have a relatively shallow gradient (Chitwood et al. 2009).

Figure 3.

One-dimensional model of adaxial–abaxial patterning generated via Hill functions and the law of mass action. Representative images of model outputs at times (A) $t = 1$ h, (B) $t = 6$ h, (C) $t = 24$ h and (D) $t = 96$ h. By $t = 96$ h, all function shapes within the model resemble experimental data for their respective gene classes.

Open in new tab Download slide

The behaviors of K, S and W also resemble those of KAN, AS and WOX seen in wild-type leaves. KAN is expressed in the two abaxial-most layers and represses transcription of AS and WOX (Nakata et al. 2012). Our model reproduces the adaxial epidermal-specific accumulation of AS, which is an outcome of epidermal transcription of AS and this direct repression by the activity of KAN (Nakata et al. 2012). WOX is transcribed in all four abaxial layers in the absence of KAN and its accumulation is detectable at approximately the P3 stage of primordial development (Nakata et al. 2012). Repression by KAN thus limits WOX accumulation to the middle two layers (Fig. 3). This is in line with experimental observations of WOX accumulation (Nakata et al. 2012). Taken together, our model recapitulates the many dynamic interactions that lead to the establishment and maintenance of the adaxial–abaxial axis [see Supporting Information—File S3] (Husbands et al. 2009, 2015).

3.2 Model interactions recapitulate biological features of small RNAs and their targets

One particularly interesting feature of small RNAs is their ability to regulate their targets via multiple mechanisms. For instance, small RNAs can clear targets from cell layers in one domain while acting as rheostats to reduce and balance (or equilibrate) their levels in the cell layers of another domain (Nikovics et al. 2006; Vaucheret et al. 2006; Skopelitis et al. 2012, 2017). We thus asked whether our model recapitulates this ‘clear and balance’ phenomenon for the HD-ZIPIII–miR166 interaction. Experimentally derived observations show the levels of targets balanced by small RNAs do not deviate by >10 % (Skopelitis et al. 2017). Therefore, let $J_{m}$ denotes the maximal interval whose endpoints are accurate to one decimal place and which yields a model containing precisely m adaxial cells and satisfying $H_{m} / H_{m - 1} \geq 0.9$ whenever $σ_{M} \in J_{m}$ ⁠. Let $I_{m}$ be defined similarly but without the condition on $H_{m} / H_{m - 1}$ ⁠. J₂ = [189.3, 257.1] as described in the Parameter Values section, and I₂ = [189.3, 266.1]. We explored model outputs for H and M for various values of $σ_{M}$ both within $I_{2}$ and outside of but near the boundary of $I_{2}$ (Fig. 4A–F). In all cases, any abaxial cell contains nearly no HD-ZIPIII mRNA. Additionally, if $σ_{M} \leq 189.2$ h⁻¹, the model results in precisely three adaxial cells, and if $σ_{M} \geq 266.2$ h⁻¹, the model results in precisely one adaxial cell. Therefore, the location of the adaxial–abaxial boundary is never ambiguous for any value of $σ_{M}$ ⁠. HD-ZIPIII mRNAs are thus cleared from abaxial cells, consistent with both the importance of the HD-ZIPIII–miR166 relationship to placement of the adaxial–abaxial boundary and the mutually exclusive nature of adaxial and abaxial cell fates (Husbands et al. 2009, 2015; Skopelitis et al. 2017).

Figure 4.

The model resolves in a biologically realistic manner across a wide range of miR166 source parameter values. (A–F) Model output of H (red) and M (green) when (A) $σ_{M} = 189$ h⁻¹, (B) $σ_{M} = 190$ h⁻¹, (C) $σ_{M} = 215$ h⁻¹, (D) $σ_{M} = 241$ h⁻¹, (E) $σ_{M} = 266$ h⁻¹ and (F) $σ_{M} = 267$ h⁻¹. The model yields precisely two adaxial cells when $190 \leq σ_{M} \leq 266$ ⁠. When $σ_{M} \leq 189$ h⁻¹ it generates three adaxial cells, and when $σ_{M} \geq 267$ h⁻¹ it generates one adaxial cell. (G) The orange curve represents $H_{3} / H_{2}$ for values of $σ_{M}$ resulting in three adaxial cells, while the blue curve shows $H_{2} / H_{1}$ for values of $σ_{M}$ resulting in two adaxial cells. Grey shaded regions correspond to $σ_{M}$ values resulting in m adaxial cells and $H_{m} / H_{m - 1} \leq 0.9$ ⁠. These unrealistic regions are short relative to the full range of $σ_{M}$ values.

Open in new tab Download slide

Next, we observed that a large range of $σ_{M}$ values result in a value of $H_{2} / H_{1}$ near 1 (Fig. 4G). Only when $σ_{M}$ is near 257.1 h⁻¹ does the value of $H_{2} / H_{1}$ start to drop noticeably below 1. Using the 10 % metric from above as a cut-off (Skopelitis et al. 2017), we can calculate $l e n g t h (J_{2}) / l e n g t h (I_{2}) = 0.8828$ ⁠. In other words, of all $σ_{M}$ values which result in the biologically expected outcome of two adaxial cells, 88.28 % also result in values of $H_{2} / H_{1}$ that are within 90 % of each other (Fig. 4G). $σ_{M}$ is thus able to fluctuate significantly within $J_{2}$ while maintaining nearly equivalent HD-ZIPIII mRNA quantities. This property also holds true at lower vales of $σ_{M}$ which produce three, rather than two, adaxial layers. Here J₃ = [127.6, 181.7], I₃ = [127.6, 189.2] and $l e n g t h (J_{3}) / l e n g t h (I_{3}) = 0.8782$ ⁠. Taken together, our model accurately recapitulates the dual regulatory mechanisms employed by small RNAs (i.e. clearing and balancing) and is robust to fluctuations in $σ_{M}$ ⁠, one of the key inputs into adaxial–abaxial patterning.

3.3 Model is robust to noise

Flat leaf production is a mechanistically challenging process, yet is remarkably robust, even in the face of genic and environmental noise. In the case of adaxial–abaxial patterning, this implies the adaxial–abaxial boundary is rarely mispositioned and that small RNA targets are cleared from one domain and balanced (or equilibrated) in the other. Accurate models of adaxial–abaxial patterning must therefore be able to buffer noise to robustly generate these outcomes.

To test the robustness of our model to noise, we introduced a new noise term to each equation in our model following standard protocols (Higham 2001). Each noise term includes a factor $a > 0$ which represents the noise amplitude; increasing a increases stochastic effects (Fig. 5). Noise implementation is discussed in greater detail in the Appendix.

Figure 5.

Noise impacts model performance. Representative images of model output at time $t = 96$ h for noise amplitude (A) $a = 0.15$ ⁠, (B) $a = 0.3$ ⁠, (C) $a = 0.45$ and (D) $a = 0.6$ ⁠. Increasing a increases stochastic effects. Only simulations with the correct numbers of adaxial and abaxial cells are shown to permit comparable visualization.

Open in new tab Download slide

a is a unitless parameter with no direct biological analog. As such, we first placed noise amplitude in a broader context by linking a to the coefficient of variation (CV) of model outputs. The CV of a set A is defined by $C V (A) = σ (A) / μ (A)$ where σ denotes the standard deviation and μ denotes the mean. Given noise amplitude a, we computed the CV for gene G in cell $C_{i}$ by defining a collection $V_{i}$ of all values of G in cell i and at time $t_{n}$ ⁠, where $t_{n}$ ranges over all discrete time points between $t = 96$ h and $t = 97$ h. The CV corresponding to cell $C_{i}$ was computed as $C V_{i} = C V (V_{i})$ ⁠. We then defined the CV for the gene G, or $C V_{G}$ ⁠, as the average value of $C V_{i}$ over all i. Finally, we defined the overall CV to be the average $C V_{G}$ over all genes G. Using these definitions, $a = 0.15$ yields a CV of 4.86 %, $a = 0.3$ yields a CV of 9.23 %, $a = 0.45$ yields a CV of 15.58 % and $a = 0.6$ yields a CV of 24.57 % (Table 2).

Table 2.

Open in new tab

Coefficient of variation values (%).

	H	S	T	M	R	K	W	Average
$a = 0.15$	3.79	7.75	3.71	3.90	5.03	4.99	4.85	4.86
$a = 0.3$	8.03	9.49	9.24	8.31	9.35	8.64	11.56	9.23
$a = 0.45$	22.06	8.29	15.65	11.95	10.07	15.24	25.81	15.58
$a = 0.6$	14.86	49.21	13.14	16.16	32.89	25.15	20.58	24.57

	H	S	T	M	R	K	W	Average
$a = 0.15$	3.79	7.75	3.71	3.90	5.03	4.99	4.85	4.86
$a = 0.3$	8.03	9.49	9.24	8.31	9.35	8.64	11.56	9.23
$a = 0.45$	22.06	8.29	15.65	11.95	10.07	15.24	25.81	15.58
$a = 0.6$	14.86	49.21	13.14	16.16	32.89	25.15	20.58	24.57

The average CV value of the model output when $96 \leq t \leq 97$ for each gene and noise amplitude parameter a. Averages in the final column are taken across all genes.

Table 2.

Open in new tab

Coefficient of variation values (%).

	H	S	T	M	R	K	W	Average
$a = 0.15$	3.79	7.75	3.71	3.90	5.03	4.99	4.85	4.86
$a = 0.3$	8.03	9.49	9.24	8.31	9.35	8.64	11.56	9.23
$a = 0.45$	22.06	8.29	15.65	11.95	10.07	15.24	25.81	15.58
$a = 0.6$	14.86	49.21	13.14	16.16	32.89	25.15	20.58	24.57

	H	S	T	M	R	K	W	Average
$a = 0.15$	3.79	7.75	3.71	3.90	5.03	4.99	4.85	4.86
$a = 0.3$	8.03	9.49	9.24	8.31	9.35	8.64	11.56	9.23
$a = 0.45$	22.06	8.29	15.65	11.95	10.07	15.24	25.81	15.58
$a = 0.6$	14.86	49.21	13.14	16.16	32.89	25.15	20.58	24.57

The average CV value of the model output when $96 \leq t \leq 97$ for each gene and noise amplitude parameter a. Averages in the final column are taken across all genes.

To measure the impact of noise amplitude on adaxial–abaxial boundary placement, we ran 1000 simulations of the model. After each simulation, we recorded the number of cells containing greater values of H than M, referred to as the number of H-adaxial cells. Likewise, we recorded the number of cells containing greater average values of T than R, referred to as the number of T-adaxial cells. The number of H-adaxial and T-adaxial cells lies between 0 and 6. For noise amplitude $a = 0.15$ ⁠, we correctly obtained two H-adaxial and two T-adaxial cells in all 1000 simulations (Fig. 6). For noise amplitude $a = 0.3$ ⁠, we correctly obtained two H-adaxial and two T-adaxial cells in 967 simulations. For noise amplitude $a = 0.45$ ⁠, we correctly obtained two H-adaxial and two T-adaxial cells in 832 simulations. Finally, for noise amplitude $a = 0.6$ ⁠, we obtained the correct adaxial–abaxial boundary placement in only 650 simulations (Fig. 6). Our model is thus able to robustly place the adaxial–abaxial boundary until a reaches a value between 0.3 and 0.45, corresponding to CVs between 9.23 % and 15.58 %.

Figure 6.

Noise affects robust positioning of the adaxial–abaxial boundary. Two-dimensional histogram showing numbers of H-adaxial and T-adaxial cells after running 1000 simulations of the model with noise amplitudes $a = 0.15$ ⁠, $a = 0.3$ ⁠, $a = 0.45$ and $a = 0.6$ ⁠. Increasing a decreases the number of simulations resulting in the correct outcome of two H-adaxial and T-adaxial cells. Note: simulations yielding H-adaxial or T-adaxial counts <1 or >3 were excluded.

Open in new tab Download slide

Next, we tested whether noise affects the ability of the model to generate cell layers with equivalent or balanced levels of transcription factors. To do this, we ran 1000 simulations of the model using various noise amplitude values. For each transcription factor G, let $A_{G} \subset {C_{1}, \dots, C_{6}}$ refers to the cells in the accumulation domain of G. Let $G_{i, j}$ refers to the value of G in cell j at time $t = 96$ h in the ith simulation. The extent to which gene G varies across cell layers in simulation i is captured by $B_{G, i} = C V ({G_{i, j} : j s a t i s f i e s C_{j} \in A_{G}})$ ⁠. Therefore, the balanced CV value for G is the average over all i of $B_{G, i}$ ⁠. For example, since H accumulates in cells $C_{1}$ and $C_{2}$ in the deterministic model, $B_{G, i} = C V ({H_{i, 1}, H_{i, 2}})$ ⁠, and we averaged these $B_{G, i}$ values over all i to compute the balanced CV value for H. These CV values are referred to as balanced CV values as they reflect the relative levels of a gene class across multiple cell layers. The null $a = 0$ noise amplitude value was included to allow comparisons with the deterministic model, and S was excluded as it accumulates in only one cell layer (Fig. 3). Consistent with expectations, increasing noise leads to higher balanced CV values for all gene classes (Table 3).

Table 3.

Open in new tab

Balanced CV values (%).

	H	R	K	W	Average
$a = 0$	0.2	1.0	0	0	0.3
$a = 0.15$	11.5	12.4	10.6	10.6	11.3
$a = 0.3$	24.8	24.9	20.7	20.7	22.8
$a = 0.45$	42.3	36.6	29.1	29.2	34.3
$a = 0.6$	52.7	48.8	38.0	39.2	44.7

	H	R	K	W	Average
$a = 0$	0.2	1.0	0	0	0.3
$a = 0.15$	11.5	12.4	10.6	10.6	11.3
$a = 0.3$	24.8	24.9	20.7	20.7	22.8
$a = 0.45$	42.3	36.6	29.1	29.2	34.3
$a = 0.6$	52.7	48.8	38.0	39.2	44.7

The average over 1000 simulations of the CV of the average mRNA values over the cells in each transcription factor expression domain. Averages are given for each protein and noise amplitude parameter a.

Table 3.

Open in new tab

Balanced CV values (%).

	H	R	K	W	Average
$a = 0$	0.2	1.0	0	0	0.3
$a = 0.15$	11.5	12.4	10.6	10.6	11.3
$a = 0.3$	24.8	24.9	20.7	20.7	22.8
$a = 0.45$	42.3	36.6	29.1	29.2	34.3
$a = 0.6$	52.7	48.8	38.0	39.2	44.7

	H	R	K	W	Average
$a = 0$	0.2	1.0	0	0	0.3
$a = 0.15$	11.5	12.4	10.6	10.6	11.3
$a = 0.3$	24.8	24.9	20.7	20.7	22.8
$a = 0.45$	42.3	36.6	29.1	29.2	34.3
$a = 0.6$	52.7	48.8	38.0	39.2	44.7

Our model is thus robust to noise with respect to placement of the adaxial–abaxial boundary, but not balancing of target levels. Specifically, adaxial–abaxial boundary placement at $a = 0.15$ is indistinguishable from the deterministic model and remains highly reproducible even at $a = 0.3$ ⁠. Only at noise amplitude values of $a = 0.45$ and $a = 0.6$ does robustness begin to break down (Fig. 6). By contrast, the robustness of target level balancing is destabilized even at low noise amplitude values. Balanced CV target levels remain relatively close to the experimentally derived 10 % CV maximum at $a = 0.15$ but are less stable than those in the deterministic model (Table 3). Robustness continues to break down as noise amplitude levels increase, and at $a = 0.6$ balanced CVs can be as high as 52.7 % (Table 3). Taken together, our results suggest the dual regulatory properties of small RNAs may have differential sensitivities to noise.

3.4 Increased substance degradation yields noise attenuation

For any substance G with Equation (8) and Solution (9), G converges to its steady state more quickly if $d_{G}$ is larger. Therefore, if noise is added to Equation (8), then increasing $d_{G}$ should buffer sudden fluctuations in G, rapidly returning the model to its steady state. For this reason, we hypothesize that degradation parameters positively correlate with robustness of the model.

This can be directly tested by comparing combinations that yield similar steady states but with different degradation parameters. We selected new parameter combinations in which all degradation parameters are either halved (×0.5) or doubled (×2). Source parameters were also halved or doubled in the respective simulations to maintain consistent levels of each gene class. The values of $σ_{M}$ and $σ_{T}$ were chosen using $J_{2}$ intervals from the deterministic model as a guide. For example, $σ_{M}$ was set at the midpoint of the interval $J_{2}$ ⁠, which correctly positions the adaxial–abaxial boundary 100 % of the time and shows no >10 % deviation of H between cells C₂ and C₁ (Fig. 4G). We ran 1000 simulations of the model with halved, standard and doubled degradation parameters for noise amplitudes of $a = 0.15$ ⁠, $a = 0.3$ ⁠, $a = 0.45$ and $a = 0.6$ ⁠, and we recorded the number of H-adaxial and T-adaxial cells for each simulation.

Standard degradation parameters with a noise amplitude of $a = 0.3$ yield the correct numbers of both H-adaxial and T-adaxial cells in 967 simulations (Fig. 7). When degradation parameters are halved, the number of correct simulations drops from 967 to 739. By contrast, doubling the degradation parameters increases the number of biologically realistic simulations from 967 to 999. Similar trends were seen for simulations run at $a = 0.15$ ⁠, $a = 0.45$ and $a = 0.6$ (Fig. 7). These data suggest degradation rate is a key factor in generating robust adaxial–abaxial patterning outcomes.

Figure 7.

Degradation parameters affect robust positioning of the adaxial–abaxial boundary. Two-dimensional histogram showing numbers of H-adaxial and T-adaxial cells after 1000 simulations of the model using noise amplitudes $a = 0.15$ ⁠, $a = 0.3$ ⁠, $a = 0.45$ and $a = 0.6$ ⁠, as well as halved, normal and doubled degradation values. Doubling degradation rates stabilized the model in the presence of increased noise, while halving degradation rates reduced robustness of model outputs. Note: simulations yielding H-adaxial or T-adaxial counts <1 or >3 were excluded.

Open in new tab Download slide

4. DISCUSSION

4.1 Inputs derived from experimental data accurately model adaxial–abaxial patterning

We have generated a robust mathematical model of adaxial–abaxial patterning. Importantly, nearly all properties of the model are estimated or informed by experimental evidence of substance interactions. One exception concerns our use of Hill functions in combination with the law of mass action (Maughan et al. 2006). We find that including these functions in small RNA-target interactions is required to generate biologically realistic outputs. Further, the Hill coefficient values required for proper model behavior are well within the range of Hill coefficients verified in other systems (Park et al. 2019). This suggests that non-linear interactions implemented using Hill functions are important components in models of adaxial–abaxial patterning.

Our model recapitulates adaxial–abaxial patterning and complements aspects of previous models (Fig. 3; see Supporting Information—File S3). For instance, a number of models have examined how the adaxial–abaxial boundary is utilized by plants to create shape (Tameshige et al. 2013; Fukushima et al. 2015; Hayakawa et al, 2016; Whitewoods et al. 2020). These models have shown how the position of the adaxial–abaxial boundary and the plane of cell division can be used to create a variety of structures. By contrast, our model describes the interactions between adaxial and abaxial determinants that position this boundary in the first place. It also suggests the interaction between HD-ZIPIII and miR166 as a key upstream determinant of robust boundary placement, and hence final organ architecture. Incorporating the substance interactions of our model might improve the predictive power of other models, providing a clearer picture of how plants develop complex morphologies.

Interactions between small RNAs and their targets have also been previously modelled. For instance, a mutually repressive relationship between small RNAs and mRNAs was shown to yield sharp on–off interfaces of target mRNAs (Levine et al. 2007). However, in this general model of small RNA-target interaction, the boundary between these two substances was allowed to occur at any point along their spatial domain. This does not accurately describe adaxial–abaxial patterning in the leaf, as adaxial and abaxial cell fates are thought to be mutually exclusive (Husbands et al. 2015). Our model improves upon this by dividing our spatial domain into discrete cell layers and restricting adaxial–abaxial boundary placement to cell walls. A mutually repressive interaction between miR166 and HD-ZIPIII gene classes has also been implemented for modelling of root patterning (Muraro et al. 2014; Miyashima et al. 2019). In these models, miR166 and HD-ZIPIII accumulate in as two inverse gradients. These predictions closely match experimental data from the root (Carlsbecker et al. 2010). Notably, when implementing the HD-ZIPIII–miR166 interaction, neither model utilized Hill functions (Muraro et al. 2014; Miyashima et al. 2019). Thus, one possible explanation for our data is that miR166 utilizes some form of cooperative behavior to regulate its HD-ZIPIII targets, but this is limited to leaf-specific contexts.

4.2 Predictions of robustness emerging from the model

The production of flat leaf architecture is a complex yet robust process. Our model predicts properties of adaxial–abaxial patterning that may contribute to this robustness. For example, by adjusting noise and degradation parameters, we found that correct resolution of the adaxial–abaxial boundary is directly proportional to the size of the degradation parameters (Fig. 7). If mRNAs involved in adaxial–abaxial patterning turn over faster than average, this could attenuate noise in the system and lead to increased robustness of outputs.

Interestingly, introduction of noise did not equally impact all aspects of the model. Positioning of the adaxial–abaxial boundary was highly robust to noise, tolerating CVs of at least 9.23 %. By contrast, even small amounts of noise destabilized the ability of the model to generate balanced levels of small RNA targets. One possible explanation is that substances and interactions required to stabilize this balancing are not included in our model. Alternatively, the dual regulatory properties of small RNA-target regulation might be under different robustness regimens. For example, the mutual exclusivity of adaxial and abaxial fates could be leveraged to increase the robustness of boundary positioning. However, this would not be a viable mechanism to ensure the robust balancing of target levels within each domain. More work is needed to assess whether these differential responses to noise reflect the biological reality.

4.3 Conclusions and future directions

Here, we have developed a model of adaxial–abaxial patterning featuring a one-dimensional fixed spatial domain. The model incorporates experimental observations of substance interaction dynamics, utilizes Hill functions in small RNA-target interactions and is highly robust to noise. Importantly, it serves as a framework to develop a more complex two-dimensional spatial model featuring growth and division of distinct cells. A two-dimensional model will enable sensitivity analyses to identify parameters that most significantly impact leaf flatness and therefore must be particularly resistant to noise. This in turn will enable us to form hypotheses of determinants and interactions that lend robustness to the mechanistically challenging process of flat leaf production. These hypotheses can then be tested using a combination of mathematical modelling and in planta molecular experiments. Our findings may also be applicable to mathematical modelling of other complex systems which generate sharp interfaces between distinct spatial regions.

SUPPORTING INFORMATION

The following additional information is available in the online version of this article—

File S1. Arabidopsis Model.m. MATLAB file containing implementation of the Arabidopsis model.

File S2. TPM Data.xlsx. Excel file containing TPM data from Cortijo et al. (2019) along with formulas to compute average TPM values per cell.

File S3. Dynamic one-dimensional model of adaxial–abaxial patterning. Movie of model outputs from $t = 0$ h to $t = 96$ h. By $t = 96$ h, all function shapes within the model resemble experimental data for their respective gene classes.

ACKNOWLEDGEMENTS

We thank Janet Best and members of the Husbands Lab for helpful discussions.

SOURCES OF FUNDING

This work was supported by National Science Foundation grant DMS-1813071, which was awarded to C.-S.C. and D. Xiu.

CONTRIBUTIONS BY THE AUTHORS

L.A. developed and implemented the model as well as wrote the manuscript. C.-S.C. provided funding and developed the model. A.Y.H. developed the model, supervised the implementation and wrote the manuscript.

CONFLICT OF INTEREST

None declared.

DATA AVAILABILITY

A MATLAB file containing implementation of the model and an Excel file containing computations used to determine values for gene TPM per cell using data from Cortijo et al. (2019) are included as Supporting Information and can be found at https://github.com/LukeAndrejek/RobustMathematicalModel.

APPENDIX

Noise Implementation

If we let $H_{j}^{i}$ denote our approximation of H in cell $C_{j}$ at time $t_{i}$ ⁠, then for $2 \leq j \leq 5$ ⁠, the finite difference formula used to compute $H_{j}^{i + 1}$ given data at time $t_{i}$ is given by

\begin{array}{l} \frac{H_{j}^{i + 1} - H_{j}^{i}}{Δ t} = \frac{D_{H}}{Δ x^{2}} (H_{j - 1}^{i} - 2 H_{j}^{i} + H_{j + 1}^{i}) + σ_{H} - d_{H} H_{j}^{i} - k_{H M} h_{H M} H_{j}^{i} M_{j}^{i} . \end{array}

We add noise to the equation for H using the formula

\begin{array}{l} \begin{array}{l} \frac{H_{j}^{i + 1} - H_{j}^{i}}{Δ t} & = \frac{D_{H}}{Δ x^{2}} (H_{j - 1}^{i} - 2 H_{j}^{i} + H_{j + 1}^{i}) + σ_{H} - d_{H} H_{j}^{i} \\ - k_{H M} h_{H M} H_{j}^{i} M_{j}^{i} + a X \sqrt{Δ t} H_{j}^{i} \end{array} \end{array}

where X is randomly sampled from a truncated normal distribution with mean 0 and standard deviation 1.

The model is implemented using an explicit numerical method since all equations in the model include non-linear terms. Such methods can struggle with stability. For this method in particular, at each time point, $Δ t$ must be small enough to satisfy a particular stability condition which depends on the values of the equations in the model at that point in time. When running the model without noise, we use a uniform value for $Δ t$ which is small enough to ensure that this stability condition is always met. When we add noise to the model, the noise term affects the stability condition. If a noise value is allowed to be chosen from a normal distribution, then given any particular value of $Δ t$ ⁠, there is always a non-zero probability that the sampled value could fail to satisfy the stability condition. However, if we randomly sample from a truncated normal distribution, then since sampled values must lie within a bounded range, we can choose $Δ t$ small enough to ensure that the stability condition will always be met. For this reason, in our simulations, X is randomly sampled from a truncated normal distribution. More specifically, we start with a normal distribution with mean 0 and standard deviation 1, and we truncate the distribution to three standard deviations.

LITERATURE CITED

Albergante

Blow

Newman

2014

Buffered qualitative stability explains the robustness and evolvability of transcriptional networks

Elife

e02863

Month:	Total Views:
December 2020	37
January 2021	107
February 2021	119
March 2021	107
April 2021	40
May 2021	57
June 2021	63
July 2021	51
August 2021	21
September 2021	43
October 2021	35
November 2021	45
December 2021	35
January 2022	43
February 2022	53
March 2022	49
April 2022	44
May 2022	32
June 2022	30
July 2022	54
August 2022	25
September 2022	41
October 2022	28
November 2022	26
December 2022	19
January 2023	31
February 2023	8
March 2023	8
April 2023	25
May 2023	20
June 2023	12
July 2023	10
August 2023	18
September 2023	38
October 2023	24
November 2023	32
December 2023	23
January 2024	28
February 2024	39
March 2024	20

Article Contents

A robust mathematical model of adaxial–abaxial patterning

Abstract

1. INTRODUCTION

2. METHODS

2.1 Identification and interaction of important substances

2.2 Equations used in model

2.3 Terms used in equations

2.4 Parameters

3. RESULTS

3.1 Adaxial–abaxial patterning can be mathematically modelled

3.2 Model interactions recapitulate biological features of small RNAs and their targets

3.3 Model is robust to noise

3.4 Increased substance degradation yields noise attenuation

4. DISCUSSION

4.1 Inputs derived from experimental data accurately model adaxial–abaxial patterning

4.2 Predictions of robustness emerging from the model

4.3 Conclusions and future directions

SUPPORTING INFORMATION

ACKNOWLEDGEMENTS

SOURCES OF FUNDING

CONTRIBUTIONS BY THE AUTHORS

CONFLICT OF INTEREST

DATA AVAILABILITY

APPENDIX

Noise Implementation

LITERATURE CITED

Supplementary data

Citations

Views

Altmetric

Email alerts

Citing articles via

Most Read

Latest

This Feature Is Available To Subscribers Only