2.2. Reaction kinetics experiments

Bovine aortic endothelial cells were cultured and transfected as previously described [

11]. We experimentally varied concentrations of bone morphogenetic protein (BMP) and matrix gamma-carboxyglutamic acid protein (MGP) presented to the cells: the BMP treatments used recombinant human BMP-4 (R&D Systems), and the MGP challenges used an N-terminally FLAG-tagged human MGP (hMGP) vector, constructed as previously described [

11], to express MGP.

The cellular response to BMP treatments (over 24 h) was recorded using luciferase assays, where the BMP-responsive luciferase reporter gene was transfected into the cells, and luciferase assays were performed as described previously and normalized to Renilla luciferase [

11].

The MGP response was measured by gene expression. Cells were plated onto 24-well plates at 3 × 10

^{4} cells per well 20–24 h prior to transfection [

12]. RNA was collected 24 h after the start of the BMP-4 treatment or transfection with the hMGP expression vector. Real-time PCR analysis was performed as previously described, and glyceraldehyde 3-phosphate dehydrogenase was used as control gene [

12,

13]. Primers and probes specific for bovine MGP were obtained from Applied Biosystems as part of TaqMan Gene Expression Assays.

2.4. Immunofluorescent staining

To visualize the expansion of cell stripes, cells were fixed in cold methanol (−20°C) for 10 min, 6 h after plating, blocked with Image-iT FX signal enhancer (Invitrogen) at room temperature for 30 min, and labelled with NMM-IIa antibody (1 : 1000, Covance, CA, USA) at room temperature for 1 h. The NMM-IIa antibodies were subsequently labelled with secondary antibodies for 30 min (Alexa Fluor 555 anti-rabbit IgG antibodies, 1 : 500; Invitrogen), and mounted by ProLong GOLD antifade with DAPI (Invitrogen). The images were acquired using a charge-coupled device (Coolsnap ES, Photometrics) equipped with an inverted microscope (Eclipse ECLIPSE TE2000, Nikon) with excitation wavelengths of 510–560 nm.

2.6. Mathematical model

Turing's original model postulated two chemical morphogens reacting and diffusing. Our experimental intervention, changing

*cell* motility, cannot be directly modelled in the original Turing model, which is a two-variable PDE in the concentrations of two chemical morphogens. Instead, the model requires a third variable to reflect cell density and the cellular processes of chemotaxis and random cell motility. Such models are well known; following the work of Keller & Segel [

15] and others [

16–

18], we modelled our system as the reaction and diffusion of a slowly diffusing chemical activator

*a*, and its rapidly diffusing inhibitor

*h*, considered as functions over a two-dimensional spatial domain (

*x*,

*y*). Our reaction terms used a version of Gierer–Meinhardt kinetics [

18]. The third variable

*N*(

*x*,

*y*) represents local cell density.

The first two equations were described previously [

9]. Briefly, in equation (2.1), the production of activator follows autocatalytic reaction kinetics and is also regulated by the inhibitor and by its own decay (at a rate

*μ*_{a}). For activation, we use a sigmoidal form

*ρ*_{a}*a*^{2}/(

*h*(1 +

*q*^{2}*a*^{2})), where

*q* is the constant for autocatalytic saturation. The inhibition of

*a* by

*h* is modelled by the

*h* term in the denominator. In equation (2.2), the production of inhibitor is proportional to

*a*^{2} and degrades at a rate

*μ*_{h}. In the first two equations, the production of activator and inhibitor is proportional to the cell density

*N.* For equation (2.3), we chose a well-known mathematical form [

7,

15] for the chemotactic term, to reflect the experimentally known fact that BMP-2 serves as a chemo-attractant agent for human vascular smooth muscle cells [

19]. Chemotactic migration following the activator gradient is regulated by a coefficient

*χ* and saturates at high levels of

*a*;

*k*_{n} is the constant for this saturation. For cell proliferation, contact inhibition of cell proliferation is a common feature by which cells restrict proliferation and cell division when the culture becomes confluent. It is particularly important for the culture of primary cells such as VMCs [

20]. In our experiment, cells reached confluence and stopped proliferating over the first week, and then accomplished the aggregation in 10–14 days. Those considerations lead to the use of logistic growth, where

*r*_{N} is the maximum rate of cell proliferation and

*N*_{c} is the cell density at confluence.

After non-dimensionalization, the equations may be expressed as:

where the dimensionless variables are

in which

*L* is the linear dimension of the domain and

*t*_{c} is the characteristic timescale of biosynthetic kinetics.

The PDEs were solved via the finite difference method. The two-dimensional domain was discretized with a uniform mesh (200 × 200). For all simulations, no-flux boundary conditions were used, and initial values of

*u* and

*v* were uniformly distributed with a 2 per cent random fluctuation, whereas the

*n* was uniformly distributed without fluctuation. In our application, the identities of the two morphogens are known: the activator

*u* is BMP-2, and the rapidly diffusing inhibitor

*v* is MGP [

9]. Parameter values were estimated from the biological literature, as described previously with some modifications [

9]. For the dynamics of morphogens, the ratio of diffusivity

*D*_{a} and

*D*_{h} was chosen as 1/200 (

*D* = 0.005) to recapture the fast diffusing inhibitor when compared with the activator, based on (i) previously measured diffusion rates of the comparable BMP homologue Decapentaplegic, (ii) previous studies of diffusivity of macromolecules in extracellular matrix, and (iii) their molecular masses [

9]. The linear degradation was conservatively estimated as 1 per cent of production rate for BMP-2 and 2 per cent for MGP, which leads to

*c* = 0.01 and

*e* = 0.02. To account for the smaller dimension of the culture plates in these experiments, we reduced the domain length,

*L*, to 1.4 cm. To allow for reduced diffusivity owing to the morphogen transportation in the newly synthesized extracellular matrix, the diffusion coefficient for the inhibitory morphogen,

*D*_{h} was estimated as 3 × 10

^{−8} cm

^{2} s

^{−1}. Using the approximate timescale of biosynthesis (1800 s), the non-dimensional scaling factor,

*γ*, becomes 35 000. Assuming cell diffusion and chemotaxis to be of the same order of magnitude, as suggested in literature [

7], we then obtained

*χ*_{0} = 0.03. Using

*r*_{N} = 0.015 h

^{−1} (based on the fact that cell number is tripled after 72 h), leads to

*r*_{n} = 322 in dimensionless form. The total time

*t** = 2 for each simulation. The constant in autocatalysis

*k* sets the saturating value at specific values of

*v*. However, as the level of

*v* (i.e. the concentration of MGP) was not known in the experiments, it is difficult to estimate the value of

*k*. As such, the suitable range of

*k* is calculated in parameter space (Turing space), in which the mechanism is driven by Turing instability. The mathematical conditions to satisfy the Turing instability are:

where the subscript denotes the first-order derivative with respect to

*u* or

*v.* The inequalities are evaluated at the steady state (

*u*_{0},

*v*_{0}). Therefore, the mathematical components in the inequalities can be written as:

The steady-state value

*u*_{0} and

*v*_{0} are functions of

*k*,

*e* and

*c*. Given

*D* = 0.005, the Turing space can be plotted as function of

*k* and

*e*/

*c* (see also [

4]). For

*e*/

*c* = 2, according to the pre-selected value of

*e* and

*c*, we can determine the range of

*k*, 0 <

*k <* 0.34, for which the parameters lie in the linearly unstable region in Turing space. In our simulation, we chose

*k* = 0.28.