Spatial variation in boundary conditions can govern selection and location of eyespots in butterfly wings

Despite being the subject of widespread study, many aspects of the development of eyespot patterns in butterﬂy wings remain poorly understood. In this work, we examine, through numerical simulations, a mathematical model for eyespot focus point formation in which a reaction-diffusion system is assumed to play the role of the patterning mechanism. In the model, changes in the boundary conditions at the veins at the proximal boundary alone are capable of determining whether or not an eyespot focus forms in a given wing cell and the eventual position of focus points within the wing cell. Furthermore, an auxiliary surface reaction-diffusion system posed along the entire proximal boundary of the wing cells is proposed as the mechanism that generates the necessary changes in the proximal boundary proﬁles. In order to illustrate the robustness of the model, we perform simulations on a curved wing geometry that is somewhat closer to a biological realistic domain than the rectangular wing cells previously considered, and we also illustrate the ability of the model to reproduce experimental results on artiﬁcial selection of eyespots.


Introduction
Eyespots, concentric bands of pigment patterning, constitute one of the most studied pattern elements on the wings of butterflies (c.f., Figure 3 for an example). Each eyespot develops around a focus, a small group of cells that sends out a morphogenetic signal that determines the synthesis of circular patterns of pigments in their surroundings. In this work, we consider a model that provides a possible mechanism underlying the determination of the number and locations of eyespots on the wing surface. The model we consider, first described in Sekimura et al. (2015), provides a mechanism that places the foci around which eyespots form in various locations on the entire wing surface. We do not address here subsequent stages of eyespot formation that occur after the development of the foci.
The model we consider is based on that of Nijhout (1990). The main novelty of the work in Sekimura et al. (2015), was to illustrate that simply changing the conditions assumed to hold at the proximal veins was sufficient to determine whether or not an eyespot formed in a given wing cell. In the present work we extend the investigations of the models proposed in Sekimura et al. (2015). We show that it is possible to determine the location of eyespots within a wing cell simply by changing the conditions that are assumed to hold at the lateral wing veins that bound the wing cell. Furthermore, we illustrate that it is possible, using a two-stage model, to recapitulate the results of artificial selection experiments in terms of selection and location of eyespots in butterfly wings.

Modelling
In this section, we describe the mathematical model for focus point formation that we consider in the present work.

Setting
As butterfly wing patterns form in two layers that are thought to be separated completely by the middle tissue (e.g., Sekimura et al., 1998), we assume that the formation of eyespots takes place in a single layer of the wing disc. Hence, we model the domain in which eyespot formation occurs as a two dimensional region. Furthermore, we assume that this two dimensional region consists of several wing cells, regions bounded by the wing veins, and we consider a region of up to seven wing cells sufficient to represent the entire surface (front or back) of the wing disc. For the sake of simplicity, we assume that each of the wing cells are of the same shape and size.
The model we consider for the formation of focus points is based on that proposed by Nijhout (1990) and consists of a reaction-diffusion system of activatorinhibitor type (Gierer & Meinhardt, 1972) posed in each wing cell with timeindependent Dirichlet boundary conditions (i.e., a source of chemicals) on the wing veins and Neumann (zero flux) boundary conditions (i.e., no flux of chemicals) at the wing margin.

Mathematical model
We denote by the number of wing cells. We denote by the 'th wing cell with boundaries , (wing margin), , , , +1 (veins) and , (proximal boundary). The boundary conditions for the activator ( 1 ) are Dirichlet (fixed) on the proximal boundary , and the wing veins , , , +1 , and Neumann (zero flux) on the wing margin , (c.f., Figure 1). The boundary conditions for the inhibitor ( 2 ) are zero flux on all four boundaries of each wing cell. The Dirichlet boundary condition on each vein , is the same for each vein. We take the initial data for both activator and inhibitor to be the positive spatially homogeneous steady state of the Gierer Meinhardt (GM) equation. Thus, our model for focus pattern formation consists of independent GM equations. The model system equations may be stated as follows: where is a diagonal matrix of positive diffusion coefficients and the reaction kinetic vector ⃗ ( ⃗) is given by 1 ( ⃗) = ( 1 1 2 2 − 2 1 ) and 2 ( ⃗) = ( 1 1 2 − 3 2 ), with 1 , 2 , 3 > 0. The choice of kinetics yields that the corresponding ODE system has a positive steady ⃗ = ( 2 / 2 , 1 3 / 2 )^. Nijhout (1990Nijhout ( ,1994 showed that the above model was capable of generating source profiles consistent with the formation of an eyespot focus within a wing cell. In Sekimura et al. (2015), we showed that changes in the Dirichlet boundary condition for 1 at the proximal boundary , alone were sufficient to determine whether or not an eyespot focus forms in a wing cell. For the proximal boundary profile, we consider two different cases firstly, prescribed boundary conditions and secondly, in order to propose a full model, we consider that the boundary profiles are themselves generated by a patterning mechanism that is posed along the entire proximal boundary, i.e., the curved surface ≔∪ , . For this one dimensional patterning mechanism, for consistency with the two dimensional model above, we consider a surface reaction-diffusion system which for illustrative purposes we choose to be the activator depleted substrate model of Schnakenberg (1979), stated as follows: is a diagonal matrix of positive diffusion coefficients, is the Laplace Beltrami operator (the analog to the usual Cartesian Laplacian on the surface) and the function ℎ ⃗⃗ (⃗⃗) is given by ℎ 1 (⃗⃗) = ( ⃗)( − 1 + 1 2 2 ) and ℎ 2 (⃗⃗) = ( ⃗)( − 1 2 2 ), with , > 0. 1 and 2 are the concentrations of two chemicals (the activator and substrate, respectively). The function can be thought of as a reaction rate and is typically taken to be constant in most studies that employ such systems to model biological pattern formation. However, if such an approach is adopted, patterns with a constant wavelength across are to be expected. In the present context, this would be insufficient to explain butterfly wing patterning in which the distribution of eyespots occurs with differing frequency in different parts of the wing. For this reason, we allow the reaction rate to be a function of space, which appears to provide sufficient freedom to generate the necessary source profiles from this one dimensional model that produce any arbitrary eyespot configuration observed on butterfly wings. The resulting model is a two-stage model for focus point formation in which the first stage corresponds to solving the Schnakenberg surface reaction-diffusion system Eq. (2.2) to steady state and in the second stage the solution 2 to this model is used to determine the proximal boundary profiles for 1 in the eyespot reaction diffusion system model Eq. (2.1) within each of the wing cells.

Computational approximation
For the approximation of the eyespot reaction diffusion system models posed within each of the wing cells, we employ an implicit-explicit finite element method developed and analysed in Lakkis et al. (2013). An advantage of such an approach is that arbitrary, potentially evolving, geometries can be considered. In particular, one does not need to assume that the wing cells are rectangular and indeed using open source meshing software, it is even possible to solve the systems on geometries obtained from image data, which may be a worthwhile extension. For the approximation of the surface reaction diffusion system, we employ the surface finite element method Dziuk and Elliott (2013). We refer to the above two references for further details on the numerical approach.

Gradients in source strength on the wing veins can determine eyespot location in the wing cell
We start by illustrating that in the eyespot focus point formation model of Section 2, it is possible to change the location of eyespots by allowing the Dirichlet boundary condition at the wing veins to vary in space. To this end, we suppose that the wing cells are trapezoidal with parallel sides corresponding to the proximal and marginal boundaries, that are chosen to be of length 1.5 and 2.5 respectively. We set the proximal boundary condition to be a convex profile of the form ( ⃗) = 2 1 (1 − sin 2 ( ( ⃗) /1.5)) where ( ⃗) is the distance from the boundary points of the proximal boundary. The boundary condition thus takes the value 2 1 at the boundary points of , and decays to 0 at the center of the proximal boundary. For the wing veins, we consider a gradient in the Dirichlet boundary condition by considering a linear boundary condition of the form ( ⃗) = 2 1 (1 − s 1 2 / 3), where 2 denotes the distance in the proximal-distal direction from the wing margin and s 1 > 0 is a parameter that governs the magnitude of the gradient. Thus the boundary condition takes the value 2 1 at the point where the vein meets the marginal boundary and decays towards the proximal boundary with slope given by s 1 > 0. The remaining parameter values we select are given in Table 1. For the discretisation we used linear finite elements on a grid with 2145 degrees of freedom (DOFs) and a timestep of 0.01. The system was solved until the discrete solution was (approximately) at steady state.

A surface reaction diffusion system model with piecewise constant reaction rate generates boundary profiles and resulting eyespot foci recapitulate those observed in artificial selection
We now report on simulations in which we illustrate that the two stage model proposed in Section 2 (see also, Sekimura et al. 2015) is capable of reproducing the differing selection of dorsal forewing eyespots observed in artificial selection experiments on Bicyclus anynana. Beldade et al. (2002) showed that, through artificial selection, it is possible to generate different phenotypes of B. anynana with either zero, one (anterior or posterior) or two forewing eyespots (anterior and posterior), c.f., Figure 3. To investigate whether our two stage model is capable of reproducing these observations, we consider a wing as shown in Figure 4. The proximal ( ) and marginal ( ) boundaries are curves corresponding to a portion of the circumference of two concentric circles of radius 9 and 12 respectively. The wing veins ( , ) are assumed to be radial and of length 3 whilst the proximal and marginal boundaries of each of the wing cells are approximately of length 1.88 and 3.35 respectively. We consider the two stage model described in Section 2. In the first stage we solve the surface reaction diffusion system with the Schnakenberg kinetics to steady state. We select Dirichlet (prescribed) boundary conditions for 1 with 1 = on one boundary and 1 = 2 at the other boundary point. For 2 we set zero-flux boundary conditions. The initial data is taken to be the steady state value for both 1 and 2 . We consider the case that the function is piecewise constant (e.g., McMillan et al., 2002), in particular we allow it to take two distinct values on either side of the midpoint (anterior-posterior) of the proximal boundary curve. The remaining parameter values we employed are shown in Table 2. After solving the Schnakenberg system to steady state we assume the Dirichlet boundary condition at the proximal boundary for the reaction diffusion system posed in each wing cell is of the form 1 ( ⃗, ) = 1.9 ̅ 2 ( ⃗) 1 ⃗ ∈ , , where ̅ 2 ( ⃗) is the spatially inhomogeneous steady state of the substrate in the Schnakenberg equation. At the veins, we set Dirichlet boundary conditions for the activator equal to twice the steady state value. The remaining parameter values are given in Table 2. We note that each wing cell in this simulation is slightly larger in area than those considered in section 4.1, and it is due to this fact that we require a slightly larger activator diffusivity, 1 , than that which was used in section 4.1.  For the numerical parameters, we used a mesh with 3927 DOFs to represent the entire wing disc. The surface reaction diffusion system was solved on the trace mesh corresponding to the boundary edges of the bulk mesh, the corresponding 1d mesh had 1793 DOFs. We used a piecewise linear finite element method for both the surface and bulk reaction diffusion systems with a timestep of 0.05 and we solved the system until the concentration profiles were (approximately) at steady state. Figure 5 shows the steady state values obtained for simulations in which we vary the value of the piecewise constant reaction rate . We see that when is zero in both the anterior and posterior, as expected the substrate concentration (that satisfies zero-flux boundary conditions) in the 1d system simply converges to a constant. Using this profile in the proximal boundary conditions for the model posed in each wing cell, we generate a wing with no foci similar to the ap case of Figure 3. If we allow to be large on one half of the proximal boundary and small on the other half, then we generate boundary profiles from the 1d system that result in a single eyespot in the half of the wing in which is large, similar to the Ap and aP phenotypes of Figure 3. Finally, if is large and constant across the entire proximal boundary, we generate a profile that leads to both the anterior and posterior foci forming as in the AP phenotype of Figure 3. The choice of Dirichlet boundary conditions for 1 leads the substrate troughs to form in the correct locations for the eventual eyespots dependent on whether they are anterior or posterior, as for zero-flux or symmetric Dirichlet boundary conditions we would expect solutions that are symmetric along the midpoint of the proximal boundary. We note that this asymmetry need not be through Dirichlet boundary conditions and could be the result of differences between individual wing cells or some other aspect which is thus far neglected in the modelling.

Discussion
In this study, we reported on further investigations of a model for the selection and distribution of eyespot foci, originally presented in the paper [Sekimura et al., (2015)]. The basic idea of the model is that whether an eyespot focus forms in a given wing cell and its eventual position in the wing cell can be determined through changing only the boundary conditions that are assumed to hold at the veins. Furthermore, we considered a two-stage model consisting of two related pattern forming mechanisms, one posed along the proximal vein and the other posed in each wing cell. The two stage model appears capable of reproducing the results of artificial selection experiments in terms of eyespot selection. A hypothesis within the two-stage model is that patterning in the first stage could be governed by a reaction-diffusion mechanism in which the reaction rate is dependent on the spatial position. Such an assumption is consistent with assuming different levels of gene activation in different regions of the wing (e.g., McMillan et al., 2002). We note however that the present model is still sensitive to changes in the parameter values and crucially, changes in the geometry. In particular, the naturally observed variations in wing cell size across butterflies appear too large for the present model to be applicable. Hence a potentially attractive avenue for future studies is to investigate Turing systems with a degree of scale invariance as has been attempted in other contexts, e.g., Othmer and Pate (1980).

Figure 5c
Steady state values of 2 and 1 for piecewise constant = 10 on one half of the wing and = 500 on the other half, corresponding to one eyespot focus on the half of the wing with increased .

Figure 5d
Steady state values of 2 and 1 for constant = 500, corresponding to two eyespot foci. Figure 5. Simulations of eyespot focus point formation using a 2-stage model. Initially a reactiondiffusion system with the Schnakenberg kinetics is solved to steady state on the curved proximal boundary using a piecewise constant value for the parameter , Dirichlet boundary conditions for