The formation of labyrinths, spots and stripe patterns in a biochemical approach to cardiovascular calcification

Calcification and mineralization are fundamental physiological processes, yet the mechanisms of calcification, in trabecular bone and in calcified lesions in atherosclerotic calcification, are unclear. Recently, it was shown in in vitro experiments that vascular-derived mesenchymal stem cells can display self-organized calcified patterns. These patterns were attributed to activator/inhibitor dynamics in the style of Turing, with bone morphogenetic protein 2 acting as an activator, and matrix GLA protein acting as an inhibitor. Motivated by this qualitative activator-inhibitor dynamics, we employ a prototype Gierer-Meinhardt model used in the context of activator-inhibitor based biological pattern formation. Through a detailed analysis in one and two spatial dimensions, we explore the pattern formation mechanisms of steady state patterns, including their dependence on initial conditions. These patterns range from localized holes to labyrinths and localized peaks, or in other words, from dense to sparse activator distributions (respectively). We believe that an understanding of the wide spectrum of activator-inhibitor patterns discussed here is prerequisite to their biochemical control. The mechanisms of pattern formation suggest therapeutic strategies applicable to bone formation in atherosclerotic lesions in arteries (where it is pathological) and to the regeneration of trabecular bone (recapitulating normal physiological development).


Introduction
Cardiovascular calcification, in atherosclerosis or valvular stenosis, is considered one of the most notorious cardiac diseases [1]. The initial atherosclerotic lesion is formed as a soft cellular and fibrous mass (atheroma) growing within the artery wall. A myocardial infarction ("heart attack") occurs when the surface of the atheroma ruptures, exposing proteins that trigger clot formation in the blood. The clots then occlude the artery, preventing blood flow to the heart muscle. Little is known about the conditions under which rupture occurs. Indeed, there is even debate about whether calcium deposits mechanically stabilize or destabilize lesions [1]. Recently it was argued that "spotty" or "speckled" patterns of calcification carry the greatest risk for plaque rupture, as opposed to uniform deposits [1]. To advance the understanding of such pathology it is therefore essential to understand the mechanism determining the patterns formed by arterial calcification.
In in vitro experiments, it was shown that cultures of vascular-derived mesenchymal stem cells differentiate and indeed spontaneously form calcified patterns out of a cellular monolayer [6]. Importantly, the morphology of the patterns was experimentally altered by applying external matrix GLA protein (MGP) (see figure 1); MGP molecules bind to active BMP-2, disabling their functionality. Thus they act as inhibitors. Assuming that the primary calcification mechanism is indeed preceded and governed by the interaction of two chemical morphogens, it is reasonable to test Turing's paradigm of morphogenesis to establish the chemical pre-patterning that shapes the pattern morphology. Alan Turing, in his seminal work on morphogenesis [7], suggested that the formation of biological patterns can be understood by means of biochemistry, that is, in the reaction-diffusion framework. In this scenario, chemicals produced by cells interact as activators or inhibitors, and diffuse through the medium at distinct rates. This can create a symmetry breaking of the uniform concentrations, a mechanism that is often called a 'diffusion-driven instability' [8]. Since then, a number of morphogens have been identified [9] and linked to pattern development [6,[10][11][12][13][14][15][16][17]. These results suggest that an understanding of the dynamics of morphogenesis can give us a relatively simple way to understand and control biological development [16,[18][19][20][21][22][23][24][25][26][27][28][29][30][31].
The Turing approach assumes a decoupling between the biochemical processes and the biomass. However, the feasibility of the Turing paradigm in physiology faces the obstacle that the primary Turing instability is linear [32,33]: the resulting pattern arises spontaneously and directly from a previously stable homogeneous condition. But the patterns that are typically observed are at large deviations from the critical conditions, and also on time scales far from the initial instability. Consequently, even if an initially homogenous biological system went through a Turing instability, it is impossible (using conventional methods) to experimentally demonstrate that fact.
In this article we discuss the role of reaction-diffusion mechanisms in promoting the mineralization patterns seen in VMSCs, as originally suggested in [6]. Using an activator-inhibitor model equation and analysis of pattern selection in the nonlinear regime, far from any initial bifurcation, Turing or otherwise, we explain previous in vitro VMSC culture experiments. Beyond that, we believe that a better description of the pattern formation phenomenon is prerequisite for biological experiments and future applications to stem cell calcification phenomena. The paper is organized as follows: in Section 2 we discuss a model equation. Next, in Section 3 we perform a one dimensional (1D) analysis to obtain the properties of periodic and localized steady states and discuss the pattern selection in the presence of multiple coexisting solutions. In Section 4, we use the latter results and secondary instabilities that operate in 2D to underline the mechanisms that lead to the formation of two dimensional patterns. We conclude and discuss the biomedical applications in Section 5.

Activator-inhibitor framework
The spatiotemporal biochemical activator-inhibitor dynamics in a VMSC culture were qualitatively modeled as follows: The monolayer of the VMSCs spontaneously expresses BMP-2 [34,35] and its inhibitor matrix GLA protein (MGP) [36,37], which diffuses more rapidly than the former [6]. With respect to local kinetics, it is assumed that BMP-2 obeys a saturated autocatalytic reaction [38] and directly promotes MGP production in a greater than linear fashion [39]. It is also assumed that both substances follow a first-order degradation. In addition, BMP-2 has a chemoattractant property [35] which is responsible for cell migration following the gradients of BMP-2. However, since: (i) the cell migration occurs at much slower time scales than the chemical diffusion, (ii) cell proliferation is relatively low [40], and (iii) neither BMP-2 nor MGP are consumed by the cells, we can neglect, to leading order, other contributions such as cell density, which indicate contributions of further, albeit minor, spatiotemporal changes of chemotactic distributions. In this approach the qualitative spatiotemporal distribution of the activator, BMP-2, emerges from interaction with the inhibitor, MGP, and thus the motion of the cells is driven by the active (MGP unbounded) BMP-2 gradients.
Following this description, we represent the active BMP-2 and MGP concentrations by a(x, y) and h(x, y) , respectively. The calcified patterns in VMSCs cultures, developed from a monolayer in a dish size of the order of centimeters, i.e., implying a large aspect ratio system (large length vs. small height) that is quasi two dimensional. An activator-inhibitor model equation of such dynamics was proposed by Gierer and Meinhardt and has often been used to study biological pattern formation [29,41] where , a h D D are diffusion constants, ρ a , ρ h are cross reaction coefficients, µ a ,µ h are degradation rates, q is the saturation constant, and A, H are source terms, respectively. We note that different variants and limits of this system have been discussed in [20,26,28,30,31,42].

Dimensionless forms, scaling, and parameters choice
For our application, we set A = 0 due to the absence of activator source in experiments [6], and rewrite (1) in a dimensionless form, by introducing dimensionless variables u = qa , v = µ a qh / ρ a , and scaling t,x by µ a −1 , D h / µ a , respectively: with S as a generalized or net inhibition source, allowed to vary.

Uniform states and bistability
Equations (2) where β ≡ 1+ 14S + S 2 , and α ≡ 108 − 72S(S + 1) + 2(S + 1) 3 . These three uniform solutions, however, exist only for a certain range of values, 0 < S < S SN ; in this range, both (0,v 0 ) and (u + ,v + ) are stable, defining a bistability region, as shown in figure 3. To leading order in S , It is important to note that it is this bistability region, and in particular the large amplitude variations in the u field, that makes possible the interspersed regions of high vs. vanishing cell densities observed in VMSC cultures [6].

Periodic solutions and localized states
In bistable systems, three primary types of stationary nonuniform solutions often arise [32,33,43]: (i) fronts connecting two uniform states (heteroclinic orbits in 1D physical space); (ii) periodic spatial patterns, also known as Turing patterns (that are limit cycles in 1D physical space); and (iii) localized states (homoclinic orbits in 1D physical space), holes or peaks superimposed on a background of a stable uniform state. The third type of state arises in systems with spatial reversibility [43,44]. In the following, we discuss the implications of periodic (ii) and localized (iii) solutions while a more general relation between solutions of type (i) and type (iii) is discussed in a companion paper [45]. First we address the temporal stability of the uniform states which is associated with non-uniform perturbations [32] ( ) , Dk Near the bifurcation onset, here Turing, wavenumber selection can be addressed using a weakly nonlinear analysis [32,46]. However, numerical integration of (2) on a large domain ( 2 / 1.38 with either periodic or reflecting boundary conditions, slightly above T S yields large amplitude periodic patterns, whose amplitudes approach 0 (0, ) v and ( , ) u v + + . This behavior indicates the presence of a subcritical bifurcation, i.e., the branch of periodic states should bifurcate first, in direction T S S < . Thus, in our case the weakly nonlinear analysis will not provide substantial insights into the nonlinear problem of pattern selection. To understand the pattern selection mechanism, we exploit first the method of spatial dynamics [44], which was found to be an effective method exploring the multiplicity of steady state solutions [43,[47][48][49][50][51][52]. To do so, we set ∂ t u = ∂ t v = 0 and rewrite (2) as a set of four first order ordinary equations where space is now treated as a time-like variable. Next, we compute the branches of steady states that bifurcate from the Turing onset which corresponds, in this specialized case (eq. (11)), to a reversible Hopf bifurcation [44], via a numerical continuation method [53]. Indeed, the instability of the ( , , , ) ( ,0, ,0) state gives rise to a subcritical bifurcation of periodic states, accompanied by a subcritical bifurcation of spatial groups consisting of a finite number of localized hole solutions [45]. Since the bifurcation is from a nontrivial state, terms 'beyond all order' in the asymptotic analysis select only two families of even parity groups, distinguished by a phase shift of π and corresponding respectively, to odd and even number of hole groups [48,54]. Notably, the simultaneous bifurcations of both periodic Turing, L = L T , and groups of localized solutions, L → ∞ , are known to arise in variational systems [47][48][49][50][51]. The bifurcation diagram of these solutions is represented in figure 3, where we use the norm (12) to distinguish among solutions, where L is the spatial period. H π and peaks P . The solutions are plotted in terms of N , see equation (12), as a function of S . All the solutions were obtained by integration of (11) via AUTO [53] on periodic domains, and their stability was determined by a numerical eigenvalue method using (2). Solid lines mark stable portions of the branch; the shaded region represents the pinning regime of the spatially bounded hole states. The 0 H and H π are branches of bounded states with a respective odd and even number of holes, each stable branch associated with a distinct number of holes which increases with N . The location of the saddle node of uniform states ( , ) u v  figure 4), at the first saddle nodes (see figure 3). The odd and even hole branches form a pinning region (see the shaded region in figure 3), where higher branches indicate increasing number of holes [45]. In variational systems, these oscillations cause an effective broadening of the Maxwell point [47][48][49][50][51] at which a heteroclinic orbit in space (Pomeau front) between the uniform and periodic solutions is time independent [55] so that the finite group of bounded holes is a connection between two heteroclinic orbits in space, i.e., a heteroclinic cycle.
The primary periodic Turing solutions . We note that there exist an infinite number of unstable wavenumbers, which for simplicity are not shown in figure 3. These solutions also bifurcate subcritically and reconnect also to the origin, but unlike the T L L = solutions, these solutions gain temporal stability after the left saddle node (see figure 3). Notably, all the periodic solutions near the left saddle nodes are an inverse image (with respect to the activator field) of the solutions (with same period) near the right saddle nodes, that is, hole-like vs. peak-like. The solutions arise first as hole-like states due to the Turing instability of the uniform nontrivial state ( ,0, ,0) u . This is demonstrated for the primary  figure 3 for details). We propose that it is this property that makes possible the chemotactic gradients that control cell migrations observed in VMSC cultures [6].
As already noted, all the periodic solution branches terminate at the onset of the transcritical bifurcation, S = 0 . The successive structure of the saddle nodes of periodic solutions above SN S S = , suggest that localized peak states, with period L → ∞ , can also be present. Indeed such homoclinic orbits in space, P , do appear but they do not bifurcate from the nontrivial branch like the periodic states. Nevertheless, by prolongation of the periodic solutions' period in the vicinity of S = 0 , we computed two distinct states, unstable (figure 4g) and stable (figure 4h), which biasymptote to 0 (0,0, ,0) v . However, unlike the possible multiplicity of the stable groups of hole solutions, only a single peak can stabilize (profile (h) in figure 4) since the double peak solution is unstable (profile (g) in figure 4), implying a repulsive interaction between two neighboring peaks. While the single type of peak has been observed in numerical integrations of (2) before [41], their origin is unclear due to the linear temporal stability of the 0 ( , ) (0, ) u v v = state (see equations (8) and (9)). To address this issue we explore the asymptotic approach of peak states to 0 (0,0, ,0) v , as x → ±∞ , which is given by where the four spatial eigenvalues are real, λ = ± E = ± 2 and λ = ± 1/ D = ± 200 . In spatially reversible systems, the formal bifurcation of localized states with monotonic tails can be either as a small amplitude (spatial) bifurcation from a uniform state [52] or by a nucleation from a global heteroclinic bifurcation, connecting in space two uniform states [45]; the latter is absent in this parameter set. In the small amplitude bifurcation case, the four real spatial eigenvalues at the bifurcation point should correspond to To qualitatively understand the peculiar connection between the periodic states and the origin of localized peaks at S = 0 , we use a projection of the four dimensional phase space onto the 2D subspace u vs. x v ∂ [44,50]. As already mentioned, in this phase space the periodic states are limit cycles; in the following we discuss only T L L = . This limit cycle persists above S = S SN , and in the absence of uniform states bistability approaches the 0 (0,0, ,0) v state as S is increased (see figure 5a). At the rightmost saddle node, the limit cycle within numerical precision crosses the trivial state 0 (0,0, ,0) v , however, the formation of a heteroclinic bifurcation between the periodic and the trivial state is impossible due to the real forms of the spatial eigenvalues [44,45,48]. This behavior, on the other hand, implies a transition to a homoclinic orbit P to the fixed point 0 (0,0, ,0) v . However, without the existence of a heteroclinic bifurcation between the two uniform states [45] a different mechanism should be considered for this special situation in which no distinction between two orbits can be made. This happens exactly at S = 0 , where uniform, periodic and peak states are simply zero; figure 5b shows the vanishing amplitude of the periodic state. Therefore the transcritical bifurcation of uniform states also serves as an effective bifurcation point for the localized peak solutions and a reason for termination of the periodic states, as demonstrated in figure 3.

Wavenumber selection in presence of multiple solutions
As has been described, the nonlinear regime exhibits the coexistence of multiple stable solutions, making it hard to foresee the wavenumber selection and the sensitivity to initial conditions; both determine the basin of attraction of the final states. To understand the mechanism, we focus here on three regions distinguished by coexisting solutions (see figure 3): (A) bounded holes and periodic states; (B) periodic states; and (C) periodic and isolated peak states.
Since in the pinning region (A), fronts between localized and periodic solutions are stationary, we refer the reader to [48], and discuss in the following only regions (B) and (C). In the linearly unstable regime T S N S S S < < , infinitesimal random perturbations around the ( , ) u v + + state grow and form periodic patterns with periodicity T L L = , due to the fastest growth of the critical Turing mode. However, other initial conditions give different results. In both regions (B) and (C), a finite amplitude spatially nonuniform initial condition that, for example, features two different length scales in two parts of the domain, results in relaxation by phase diffusion [56,57]: the asymptotic periodicity is a compromise value between two initial periodic states, implying dispersion of the front separating two domains with different periodicities rather than invasion or phase slips, see figure 6. When we initiated a domain with S SN two different stable solutions coexisting in regions (B) and (C), it then evolved into a single final pattern, whose length scales were intermediate, which are also stable solutions of (2) (figures 6a and 6b, respectively).

Mechanisms of pattern formation in 2D
It is well-known that reaction diffusion systems can exhibit transverse and curvature-induced instabilities in 2D [32,33,58], so that labyrinthine or spotted patterns can form respectively, through secondary zigzag or varicose instabilities of stripes [59] or through a transverse instability of axisymmetric spots [60]. Thus, we have performed an extended numerical investigation to obtain the secondary pattern formation mechanisms of the diverse patterns. First, we have numerically calculated the regions of zigzag, varicose and curvature instabilities of periodic stripes, isolated stripes, and isolated spots, respectively. The distinct instability types and their ranges indicated in figure 7 and respectively shown in figure 8. We have exploited the 1D periodic and localized solutions obtained above to construct the respective periodic (figure 8a) and isolated stripe (figure 8a) solutions in 2D while the spot (figure 8c) solutions were obtained by a relaxation from an axisymmetric construction based on the initial localized peak or hole state size. We do not show the transient evolutions from spots toward labyrinthine patterns in the regime T S N S S S < < , and refer the reader to discussions of similar phenomena in the context of isolated stripes [28] and isolated spots [60].  One of the primary interests in the experimental context is the control of the chemical pre-pattern under various concentrations of the inhibitor source. Thus, to demonstrate the effect of these instabilities on the asymptotic pattern formation, we have at first integrated equation (2) in 2D above the Turing onset, starting from infinitesimal random perturbations around the uniform (u + ,v + ) state; as expected due to the zigzag instability (see figure 8a), a labyrinthine pattern was formed (see inset (a) in figure 7). Next, we used this pattern as an initial state for other values of S , the results are shown in insets (b-g) in figure 7. The domain of labyrinthine-type patterns ranges from the bistability onset at S = S SN (inset (b) in figure   7) to the 2D pinning region (inset (c) in figure 7), however only for S T < S < S SN do they form spontaneously, due to the Turing instability of the uniform (u + ,v + ) state. In a small region above S SN , single spots of peak type are still unstable (see figure 8e), so that straight stripe fragments may form. In case of initial labyrinthine pattern (inset (a) in figure 7) this leads to the formation of periodic states consisting of both spots and stripe fragments (inset (d) in figure 7). As a single spot stabilizes (see figure 8f), the pattern becomes ultimately spotted. However, according to figure 3, the pattern may admit distinct periodicities under other initial conditions; direct numerical 2D computations support this prediction, as insets (e-f) in figure 7 show. We note that both states (e-f) are not asymptotic and approach asymptotically to a hexagonal symmetry, however, since the interaction between localized states is weakly repulsive, this process is very slow. For small S values (that is, in the 2D pinning region [61,62]), we obtain an image inverse to peak spots, that is, distributed hole spots (inset (g) in figure 7). In the absence of repulsive interactions between neighboring holes, see insets (h-i) in figure 7, the hole spots support coexisting clusters, as shown in inset (g) in figure 7.

Conclusions
In this paper, we theoretically studied the mechanisms of pattern formation that give rise to many steady state patterns in a biologically oriented activator-inhibitor model [29,41]. We chose a version with an inhibition source term, and a saturated autocatalytic local kinetics. These resulted in a bistability regime of uniform solutions, distinguished by finite and zero activator concentrations (see equation (6)). While different forms of the Gierer-Meinhardt model have been studied, and a number of interesting qualitative phenomena have been found [20,26,28,30,42], we present here a generalized and detailed view of the multiple periodic and localized states that occur in this model, as well as their basins of attraction. In particular, we show that the majority of patterns for high inhibitor source values form subcritically, that is, by finite amplitude perturbations above the trivial state 0 (0, ) v , as shown in figures 3 and 7. First, we exploited the powerful tool of spatial dynamics coupled with numerical continuation, and incorporated with temporal stability, to identify the primary instabilities and the possible stable and unstable 1D solution branches as a function of the inhibition source, S (see figure 3). We used this analysis to calculate the secondary instabilities and the parameter regions in which 2D patterns, such as labyrinths, mixtures of spots and stripe fragments, and bounded vs. isolated spots, can form (see figure 7). The knowledge of the nonlinear pattern selection mechanisms far from the Turing onset is therefore an important missing link in the context of robustness and pattern control, which up to now was considered as one of the weak points of the activator-inhibitor approach in biological pattern formation [63]. The qualitative form of the model, and in particular the parameter values, were made to keep fidelity with our main application, which is the formation of calcified patterns by vascular mesenchymal stem cells (VMSCs) [6]. These are the cells that are thought to form bone in atherosclerotic tissue, the primary process that has been linked to atherosclerotic calcification [2]. The approach we employ is associated with the primary symmetry breaking mechanism of a cellular monolayer. This assumption is justified by the fact that cells do not consume the morphogens, allowing a qualitative decoupling between the biochemistry and the biomass. Thus, we have excluded from this model any higher order contributions, such as the contributions of changing cell densities, which may support chemotaxis-based mechanisms [30]. Within these lines, cell density is expected to have a minor role in the initial qualitative pattern selection, and only quantitatively modify the pattern architecture at later time stages. The extension of the Gierer-Meinhardt model to include the contributions of cell density is beyond the scope of this paper, and will be discussed elsewhere.
Importantly, a spontaneous formation of labyrinthine and spotted patterns (with roughly hexagonal symmetry) has indeed been observed in in vitro experiments in VMSC preparations, as shown in figure 1, and discussed in more detail in [6]. In these experiments the inhibition source was altered by external addition of MGP, the inhibitor of the activator BMP-2. As in our analysis, labyrinthine patterns have been found at lower concentrations of added MGP (corresponding to lower S values) while at higher MGP concentrations (corresponding to higher S values) periodic spotted patterns were developed. The good qualitative agreement between the theoretical analysis and the experimental observations suggests a novel strategy for the biochemical control of calcified patterns, via the framework of activatorinhibitor dynamics.
There are several potential applications of this strategy. First, as noted above, the spotted patterns of calcification seen in culture dishes are considered a good model for the formation of spotty atherosclerotic calcification in humans [1]. If these calcified deposits become confluent, it has been conjectured that the resulting solid mass is expected to mechanically stabilize the adjacent lesion, while if the calcification remains 'spotty', there is an increased risk of rupture and myocardial infarction [1]. We believe that the theoretical prediction of patterns will shed light on design of future experiments and allow also identification of the proper biochemical conditions for spot nucleation phenomenon (peak states), which according to the theory can not be observed spontaneously and require specific initial conditions. A second application is to the formation of bone tissue (by the same types of cells). Trabecular bone has a dense labyrinthine and hole architecture, presumably formed by the same morphogenetic processes, so our analysis may also be relevant to biochemical based methods of bone regeneration.