Turing mechanism underlying a branching model for lung morphogenesis

The mammalian lung develops through branching morphogenesis. Two primary forms of branching, which occur in order, in the lung have been identified: tip bifurcation and side branching. However, the mechanisms of lung branching morphogenesis remain to be explored. In our previous study, a biological mechanism was presented for lung branching pattern formation through a branching model. Here, we provide a mathematical mechanism underlying the branching patterns. By decoupling the branching model, we demonstrated the existence of Turing instability. We performed Turing instability analysis to reveal the mathematical mechanism of the branching patterns. Our simulation results show that the Turing patterns underlying the branching patterns are spot patterns that exhibit high local morphogen concentration. The high local morphogen concentration induces the growth of branching. Furthermore, we found that the sparse spot patterns underlie the tip bifurcation patterns, while the dense spot patterns underlies the side branching patterns. The dispersion relation analysis shows that the Turing wavelength affects the branching structure. As the wavelength decreases, the spot patterns change from sparse to dense, the rate of tip bifurcation decreases and side branching eventually occurs instead. In the process of transformation, there may exists hybrid branching that mixes tip bifurcation and side branching. Since experimental studies have reported that branching mode switching from side branching to tip bifurcation in the lung is under genetic control, our simulation results suggest that genes control the switch of the branching mode by regulating the Turing wavelength. Our results provide a novel insight into and understanding of the formation of branching patterns in the lung and other biological systems.


Introduction
The Mammalian lung is a striking example of organs that develop through branching morphogenesis. During lung morphogenesis, two primary forms of branching, side branching and tip bifurcation, which occur in sequence, have been identified [1]. The switch of branching mode from side branching to tip bifurcation is postulated to be under genetic control [1,2].
To investigate how genes work to generate these patterns, a mathematical model [3] derived from the Gierer-Meinhardt activator-inhibitor model [4] was used in our previous a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 study [5]. We demonstrated a mechanism through which the interaction of biological morphogens creates branched structures in the lung. The cascades of branching forms that have been observed in the lung, including side branching and tip bifurcation, were successfully reproduced by the branching model. Although the biochemical mechanism-the interaction of morphogens-provides an elegant explanation of lung branching morphogenesis, the mathematical mechanism underlying the branching patterns needs to be further investigated. For example, the branching mode switch between side branching and tip bifurcation can be controlled by a key parameter related to consumption by cells in the simulation of the model; however, it is not easily explained by the interaction of morphogens. Mathematical studies focus on the dynamical behaviors of mathematical models [6][7][8][9]. However, there is lack of bridge between branching morphogenesis and mathematical mechanism. Based on the branching model, we investigate the mathematical mechanism underlying lung branching pattern formation in this paper.
In our previous study of the dynamics of side branching and tip bifurcation [10], we showed that Turing instability occurs in the branching patterns. Turing instability can induce spatial patterns in the models, such as spots, stripes, hole patterns, and more complicated patterns, which is applied to modeling biological patterning phenomena in fish skin, terrestrial vegetation, sea shells, and others [11][12][13][14]. To reveal the mathematical mechanisms underlying branching patterns, we conducted Turing instability analysis.
In this paper, we decoupled an activator-inhibitor model from the branching model and performed simulations of the two models to obtain Turing patterns and branching patterns. Our simulation results show that Turing instability occurs at the growing tip of the branching patterns. The Turing patterns underlying the branching patterns are spot patterns. The spot patterns are in the form of concentration peaks, leading to branching patterns, with a local activator concentration peak formed and moving ahead of the growing tips. This indicates that the local morphogen concentration peak plays a key role in the growth of branching. Furthermore, a sparse spot pattern underlies the tip bifurcation patterns, while a dense spot pattern underlies the side branching patterns. The dispersion relation analysis shows that the wavelength of the spot patterns acts on the branching structures. As the wavelength decreases, the spot patterns change from sparse to dense, the rate of tip bifurcation decreases and side branching eventually occurs instead. A sufficient wavelength is required for the occurrence of tip bifurcation, while an insufficient wavelength provides favorable conditions for side branching. Our results suggest that genes control the branching structures in the lung by regulating the Turing wavelength. Our results provide a fresh insight into and understanding of the formation of branching patterns in the lung and other biological branching systems.

Methods
The branching model we used in this paper is defined by Eqs (1)-(4). The four variables in the model equations are: activator A, inhibitor H, substrate S, and cell differentiation state Y.
In the branching model, the term cA 2 S H in the first equation represents that activator A is upregulated by itself in autocatalytic reaction kinetics at rate c, with the dependence on substrate S, and inhibited by inhibitor H. cA 2 S in the second equation represents the catalyzed effect of H by A. ρ A Y and ρ H Y represent A and H are secreted by differentiated cells Y at rates ρ A and ρ H . c 0 in the third equation represents S is produced at rate c 0 . −εYS represents S is consumed by differentiated cells Y at a rate ε. dA þ Y 2 1þfY 2 represents a high A-concentration induces irreversible cell differentiation (the Y-concentration goes from low to high). −μA, −vH, −γS, and −eY represent A, H, S, and Y decay in a first-order reaction at rates μ, v, γ, and e. A, H, and S are assumed to be diffusible, with diffusion coefficients D A , D H , and D S , respectively.
This branching model is derived by an activator-inhibitor system (Eqs 1 and 2) with dependence on the substrate. Thus, the branching pattern formation described in the model is divided into two processes: the formation of the local pattern on the stalk and the spatial extension of the stalk. The former is generated by the activator-inhibitor system, which exhibits Turing instability, and the latter results from the dependence of the activator-inhibitor system on the substrate. To explore the Turing instability underlying the branching model, a scheme was performed according to the following steps, as depicted in Fig 1. Step 1. Decoupling. Decouple the activator-inhibitor model from the branching model, with S and Y as the parameters. We used the decoupling method described in the literature [10].
Step 2. Calculating a crescent-shaped Turing region. Calculate the S-Y parameter space of the activator-inhibitor model for the Turing instability (see Appendix for the Turing instability analysis).
Step 3. Plotting the differentiation trajectory of a cell. Extract (S, Y) pairs of a cell in the branching system with cell differentiation, then plot an SY-curve for the differentiation trajectory of the cell.
Step 4. Turing state selection and simulation. Select a point on the trajectory within the Turing region as the values of parameters S and Y and perform simulation of the activator-inhibitor model to obtain the Turing-type pattern underlying the branching pattern.

Numerical simulation
We performed numerical simulations of both the branching model and the activator-inhibitor model to obtain the branching patterns and the related Turing patterns. For the branching model, we set the values of the parameters according to the literature [5]. For the activatorinhibitor model, simulations were performed on a 200×200 grid with periodic boundary conditions, and the parameter values and initial values of the variables were set according to the branching system. Starting from a randomly perturbed uniform initial condition, the simulation of the activator-inhibitor model was stopped when the stationary spatial structure was formed.

Turing spot patterns underlying the branching patterns
To explore the Turing patterns underlying the branching patterns, we calculated the S-Y parameter space for the Turing instability to obtain a crescent-shaped Turing region and recorded the cell differentiation trajectories of both the tip bifurcation and side branching patterns. We then performed simulations of their underlying Turing patterns. The simulation results are shown in Fig 2. For the tip bifurcation patterns (Fig 2Aa), the underlying Turing patterns are spot distributions (Fig 2Ab). The spot patterns are at points on the cell differentiation trajectories of the tip bifurcation patterns within the Turing region (Fig 2Ac).The points refer to a cell state where cells are located at the growing tips of the tip bifurcation patterns. This means the Turing instability occurs at the growing tips of the tip bifurcation patterns.
To further investigate the effect of the spot patterns on the tip bifurcation patterns, we explored their activator distribution since the activator and inhibitor interaction plays a key role in pattern formation. In Fig 2C, the tip bifurcation patterns always exist with a local activator concentration peak that is formed and moves ahead of the growing tips (Fig 2Ca and 2Cb). The spot patterns are in the form of peaks of the activator concentration (Fig 2Cc). This indicates that the Turing instability affects the branching growth. The structure of the spot pattern leads to the activator concentration peaks formed at the growing tips of the tip bifurcation patterns.
We then addressed how the Turing spot patterns underlie the tip bifurcation patterns. To interpret the mechanism, we explored the morphogen concentration in the Turing patterns because it plays an important role in cell growth [15,16]. There are other typical structures of Turing patterns in addition to spots, such as stripes and holes. Fig 3 shows the activator concentration of those patterns. Spot patterns exhibit local high concentration peaks, while stripe and hole patterns are shown in a gentle concentration distribution. The spot patterns have a much higher activator concentration gradient than the stripe and hole patterns. This indicates that a high local morphogen concentration is required for tip bifurcation growth. The high concentration peaks at the growing tips of the tip bifurcation patterns caused by the Turing instability stimulate the outward extension of the tips and outward growth of the branches.
With respect to the side branching patterns (Fig 2Ba), the underlying Turing patterns also have spot distributions (Fig 2Bb). The spot patterns are at points on the cell differentiation trajectories of the side branching patterns within the Turing region (Fig 2Bc). The points refer to the cell state where cells are located at the growing tips of the side branching patterns, which means Turing instability also occurs at the growing tips for the side branching patterns. Fig 2D shows that a local activator concentration peak is formed and moves ahead of the growing tips of the side branching patterns (Fig 2Da and 2Db), which is consistent with the spot patterns in  the form of peaks of the activator concentration (Fig 2Dc). Those results are the same as the case of the tip bifurcation patterns. However, the spot patterns underlying the side branching patterns are much denser than those corresponding to the tip bifurcation patterns.

Spot density of the Turing spot pattern varies for branching structures
We observed that the structure of the Turing patterns underlying the branching patterns is spots. Furthermore, the spot density of the Turing patterns varies for the tip bifurcation and side branching patterns. Next, we investigated the connection between the spot density of the Turing spot patterns and the branching structures. In the simulation, both branching structures, tip bifurcation and side branching, can be generated by modifying a single parameter, ε (in Eq (3), the consumption rate of substrate by Y cells). In this way, we obtained the branching patterns and the underlying spot patterns, as shown in Figs 4 and 5.
Sparse Turing spot patterns underlying the tip bifurcation patterns. We set the tip bifurcation pattern shown in Fig 2Aa for a given ε. For convenience of comparison, we show the pattern in Fig 4B and its underlying spot pattern in Fig 4E. We then increased ε and obtained a tip bifurcation pattern with an increasing bifurcation rate, and the underlying spot pattern with a decreasing number of spots was observed (Fig 4A and 4D). Subsequently, we decreased ε, and another tip bifurcation pattern with a decreasing bifurcation rate were obtained, while the underlying spot pattern with an increasing number of spots was observed (Fig 4C and 4F).
For the tip bifurcation patterns, the underlying Turing patterns are sparse spot patterns. Tip bifurcation occurs at a decreasing rate with increasing number of spots in the spot patterns as ε decreases.
Dense Turing spot patterns underlying the side branching patterns. When ε is below a certain value, side branching patterns emerge rather than tip bifurcation patterns. We set the side branching pattern shown in Fig 2Ba for a given ε; we show the pattern in Fig 5B and its underlying spot pattern in Fig 5E. We then increased ε within the range for side branching patterns, and a side branching pattern was obtained with slightly increasing spatial interval between branches, while the underlying spot pattern was observed with a slightly decreasing number of spots (Fig 5A and 5D). Subsequently, we decreased ε, and another side branching pattern was obtained with a slightly decreasing spatial interval between branches and more outward growing branches, while the underlying spot pattern was observed with a the slightly increasing number of spots (Fig 5C and 5F).
For the side branching patterns, the underlying Turing patterns are dense spot patterns. As ε decreases, more side branches are produced, the spatial interval between branches decreases, and the number of spots in the underlying spot patterns increases.

Turing wavelength regulates the branching structures
Turing patterns are characterized by a critical wavelength [17]. To elucidate the phenomenon of distinct spot densities of the spot patterns underlying the tip bifurcation patterns and side branching patterns, we further explored the critical wavelength of the spot patterns by dispersion relation analysis. The dispersion relations shown in Fig 6 describe a function of Re(λ) that depends on wavenumber k, where λ is the eigenvalue with the largest real part (see Appendix for obtaining the dispersion relations). The wavelength is calculated by dividing 2π by the critical wavenumber at which the maximum value of λ occurs.
In Fig 6, we show the dispersion relations for the spot patterns underlying the branching patterns depicted in Figs 4 and 5 and present a comparison of the wavelength (2π/wavenumber) in sequence. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot patterns increases (Fig 4), the dispersion relations illustrate that the wavelength decreases (Fig 6A, green curves; Fig 6B, green bars). When the branching mode switches from tip bifurcation to side branching with the underlying spot pattern change from a sparse to dense distribution, the wavelength decreases further (Fig 6A,  green to orange curves; Fig 6B, green to orange bars). For the side branching patterns, when the spatial interval between branches decreases with more outward growing branches, the number of spots in the underlying dense spot pattern increases slightly, and the wavelength decreases slowly (Fig 6A, orange curves; Fig 6B, orange bars). These data suggest that as the wavelength decreases, the number of spots in the spot patterns increases, and tip bifurcation occurs at a decreasing rate. When the wavelength decreases below a certain value, side branching occurs rather than tip bifurcation.
To investigate the effect of the wavelength on the branching patterns for different Turing regions, we varied parameter ρ H (in Eq (2), inhibitor secreted by cells) to explore the different Turing regions, since ρ H is a key factor in the branching patterns in the simulation [5] and is a parameter in the activator-inhibitor model (Eqs (1) and (2)). As we have already analyzed the Turing region for ρ H = 0.0001, four Turing regions for two smaller and two larger ρ H values are selected for the analysis. Through dispersion relation analysis (Fig 8C), we show the wavelength (Fig 8D). Trends similar to those for the Turing region with ρ H = 0.0001 were observed. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot patterns increases (Fig 8B), the wavelength decreases (Fig 8D, green bars). When the branching mode switches from tip bifurcation to side branching with the underlying spot patterns changing from a sparse to a dense distribution (Fig 8B), the wavelength decreases further (Fig 8D, green to orange bars). For the side branching patterns, when the spatial interval between branches decreases and more branches grow outward, the number of spots in the underlying dense spot patterns increases slightly (Fig 8B), and the wavelength decreases slowly (Fig 8D, orange bars). In addition to the similar trends, an interesting hybrid branching pattern (Fig 8B4), which mixes tip bifurcation and side branching, is generated. Both the number of spots and the wavelength of the underlying spot patterns are between those corresponding to the tip bifurcation and side branching patterns (Fig 8B4 and 8D, blue bars).
For the Turing region for ρ H = 0.00007, similar trends are observed in Fig 9. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot patterns increases (Fig 9B), the wavelength decreases (Fig 9D, green  bars). When the branching mode switches to side branching with the underlying dense spot patterns (Fig 9B), the wavelength decreases further (Fig 9D, green to orange bars). For the side branching patterns, when the branches grow closely and the number of spots in the underlying dense spot patterns increases slightly (Fig 9B), the wavelength decreases slowly (Fig 9D, orange bars).
Similarly, a hybrid branching pattern also emerges, and both the number of spots and the wavelength of the underlying spot patterns are between those corresponding to the tip Turing mechanism underlying a lung branching model  bifurcation and side branching patterns (Fig 9B4 and 9D, blue bars). However, Fig 9B4 shows that it is not easy for side branch to emerge from the tip branching structure.
Similar trends were found in the Turing region for ρ H = 0.00013 in Fig 10. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot patterns increases (Fig 10B), the wavelength decreases (Fig 10D,  green bars). When the branching mode switches to side branching with a dense underlying spot pattern (Fig 10B), there is an evident decrease in the wavelength (Fig 10D, green to orange bars). For side branching patterns, when branches grow close and the number of spots in the underlying dense spot patterns increases slightly (Fig 10B), the wavelength decreases slowly (Fig 10D, orange bars).
However, no hybrid branching patterns were observed in the Turing region for 0.00013. Similar trends were observed in Fig 11 in the Turing region for ρ H = 0.00015. When the bifurcation rate in the tip bifurcation patterns decreases and the number of spots in the underlying sparse spot patterns increases (Fig 11B), the wavelength decreases (Fig 11D,  green bars). When the branching mode switches to side branching with an underlying dense spot pattern (Fig 11B), the wavelength decreases greatly (Fig 11D, green to orange bars). In the side branching patterns, when branches grow close and the number of spots in the underlying dense spot patterns increases slightly (Fig 11B), the wavelength decreases slowly (Fig 11D, orange bars).
Additionally, there are no hybrid branching patterns observed in the Turing region for ρ H = 0.00015.
The simulation results demonstrate that the effects of the wavelength on the branching patterns have similar trends in different Turing regions. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot pattern increases, the wavelength decreases. When the branching mode switches from tip bifurcation to side branching and the underlying spot pattern changes from a sparse to dense distribution, the wavelength decreases further. For the side branching patterns, when the spatial interval between branches decreases and the number of spots in the underlying dense spot pattern increases slightly, the wavelength decreases slowly.

Discussion
Our simulation results demonstrate that a local high morphogen concentration and the Turing wavelength play important roles in pattern formation in the branching model.
In the branching patterns, the growing tips exhibit Turing instability. We show that it is the Turing spot patterns underlying the branching patterns. The spot patterns are in the form of concentration peaks, which results in a local morphogen concentration peak formed at the tips in the branching patterns. The local morphogen concentration peak is unstable and induces tip expansion into the free space, causing branches to grow. This result is in agreement with the in vitro experimental results of Hagiwara et al [18], who showed that a high cell concentration gradient is required for cell branching in the lung.
Furthermore, we found that the spot density of the spot pattern varies for branching structures. A sparse spot pattern underlies the tip bifurcation patterns, while a dense spot pattern underlies the side branching patterns.
The dispersion relation analysis shows that the wavelength of the spot patterns affects the occurrence of tip bifurcation and side branching. For the tip bifurcation patterns, when the bifurcation rate decreases and the number of spots in the underlying sparse spot pattern increases, the wavelength decreases. When the branching mode switches from tip bifurcation to side branching and the underlying spot pattern changes from a sparse to a dense distribution, the wavelength decreases further. For the side branching patterns, when the spatial interval between branches decreases as more branches grow, the number of spots in the underlying dense spot pattern increases slightly and the wavelength decreases slowly.
The simulation results suggest that when the wavelength decreases and the number of spots in the spot patterns increases, tip bifurcation occurs at a decreasing rate. When the wavelength decreases below a certain value, the spot patterns shift to a dense distribution, no tip bifurcation occurs and side branching is observed. An insufficient wavelength impedes tip bifurcation but provides favorable conditions for side branching.
Branching patterns and the Turing patterns are two types of patterns in mathematical biology. Our work contributes to correlating the formation of branching patterns with Turing patterns. Although we demonstrate the connection between spot patterns and branching patterns, little is known about other Turing patterns, such as stripe patterns and hole patterns. The dispersion relation analysis shows that the wavelength affects the branching pattern, and the trend of how the wavelength affects the branching structures is revealed. However, the exact mechanism remains to be explored.
Nevertheless, our work reveals the Turing mechanism underlying the branching patterns. In our previous study [5], we demonstrated that the branching mode can be changed from tip bifurcation to side branching by varying the parameter ε. Our results in this paper further show that ε controls the branching mode switch because it regulates the Turing wavelength. In the experimental work, the branching mode changes during lung development are shown to be controlled by genes. Our results further suggest that the branching mode switch in the lung is a result of genes regulating the Turing wavelength, similar to a previous study [19], which found that gene modulation of digit patterning involves a Turing mechanism. Our work provides a fresh insight into and understanding of the formation of branching patterns in the lung and other biological branching systems.