A Model for Selection of Eyespots on Butterfly Wings

Unsolved Problem The development of eyespots on the wing surface of butterflies of the family Nympalidae is one of the most studied examples of biological pattern formation.However, little is known about the mechanism that determines the number and precise locations of eyespots on the wing. Eyespots develop around signaling centers, called foci, that are located equidistant from wing veins along the midline of a wing cell (an area bounded by veins). A fundamental question that remains unsolved is, why a certain wing cell develops an eyespot, while other wing cells do not. Key Idea and Model We illustrate that the key to understanding focus point selection may be in the venation system of the wing disc. Our main hypothesis is that changes in morphogen concentration along the proximal boundary veins of wing cells govern focus point selection. Based on previous studies, we focus on a spatially two-dimensional reaction-diffusion system model posed in the interior of each wing cell that describes the formation of focus points. Using finite element based numerical simulations, we demonstrate that variation in the proximal boundary condition is sufficient to robustly select whether an eyespot focus point forms in otherwise identical wing cells. We also illustrate that this behavior is robust to small perturbations in the parameters and geometry and moderate levels of noise. Hence, we suggest that an anterior-posterior pattern of morphogen concentration along the proximal vein may be the main determinant of the distribution of focus points on the wing surface. In order to complete our model, we propose a two stage reaction-diffusion system model, in which an one-dimensional surface reaction-diffusion system, posed on the proximal vein, generates the morphogen concentrations that act as non-homogeneous Dirichlet (i.e., fixed) boundary conditions for the two-dimensional reaction-diffusion model posed in the wing cells. The two-stage model appears capable of generating focus point distributions observed in nature. Result We therefore conclude that changes in the proximal boundary conditions are sufficient to explain the empirically observed distribution of eyespot focus points on the entire wing surface. The model predicts, subject to experimental verification, that the source strength of the activator at the proximal boundary should be lower in wing cells in which focus points form than in those that lack focus points. The model suggests that the number and locations of eyespot foci on the wing disc could be largely controlled by two kinds of gradients along two different directions, that is, the first one is the gradient in spatially varying parameters such as the reaction rate along the anterior-posterior direction on the proximal boundary of the wing cells, and the second one is the gradient in source values of the activator along the veins in the proximal-distal direction of the wing cell.


Introduction
Butterfly wing color patterns are among the most spectacular and remarkable examples of patterning in biology. For more than a century, they have attracted much attention from experimentalists and theoreticians alike. One of the most studied color patterns on butterfly wings is the eyespot (Fig 1) that may play a central role in interactions with predators. The formation of eyespots has been the subject of studies in molecular and developmental genetics (e.g., [1,2,3]), evolution, physiology (e.g., [4,5]), ecology (e.g., [6,7,8]), and theoretical biology (e.g., [9,10,11]). These studies, however, have focused on the formation mechanism of a single eyespot located at a specific position on the wing surface. Several species of butterflies, however, have many eyespots on their wing surface. The number, size, shape, pigmentation and precise position of these eyespots are extremely diverse and are typically species-specific. In order to fully understand the evolution and diversity of eyespot patterns, it is necessary to analyze the mechanism that governs the formation of these different pattern elements.
In this paper, we focus on the mechanism underlying the determination of the number and locations of eyespots on the wing surface. 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. Our paper is concerned with the mechanism that places these foci in various locations on the wing surface. This fundamental process constitutes the first of three developmental steps of eyespot formation (see Section 2 for details). We do not consider the mechanism behind the determination of the size, shape and pigment patterns around the foci, which occurs later in the developmental stages. The number and locations of foci would have undergone considerable evolution during the diversification of the butterflies. Our objective in this article is to propose a model that determines the global distribution of foci in the overall venation system of the wing.

General features of eyespot formation
The formation of wing color patterns including eyespot patterns is a spatially two-dimensional phenomenon that takes place in the single layer of cells that makes up each surface of the wing [e.g., 12]. The butterfly wing begins its development as a wing imaginal disc in the larva. The wing imaginal disc is transparent and colorless throughout the larval and early pupal stages of development. Antibody and mRNA fluorescence techniques for the several developmental genes have revealed existence of a developmental pre-pattern on the wing disc, which predicts the color pattern of the adult butterfly wing (e.g., [3]).
For the specific case of eyespot formation, the formation mechanism is thought to consist of the following three developmental stages (e.g., [13, 14]) i. The first stage is the determination of the location of the signaling center, i.e., "the eyespot focal cells", from which some signaling chemicals, or morphogens originate.
ii. The second stage is the spreading out of morphogens into the surroundings of the focus cells through diffusion and activation of corresponding genes (e.g., Dll, engrailed), which establish a pattern of concentric rings of gene expression that constitutes the pre-pattern for pigment synthesis.
iii. The third stage is the activation of the pigmentation genes (e.g., DDC, GTP-CH1, cinebar) that cause the synthesis of species-specific pigments as a set of concentric colored rings we recognize as an eyespot. The focus cells are typically pigmented white and form the "pupil" of the eyespot on the adult wing.
The target problem in this paper pertains to the first developmental stage described above: the formation and positioning of the foci. We will not consider the growth of the wing disc as part of the modeling as this is assumed to occur on a longer timescale influencing only the second and the third stages.

A mechanism for selection of eyespot focus points in the wing disc
Although there is experimental data on the development of eyespot foci, little is known about the mechanism that determines the number and locations of focus points in the entire wing disc. As seen in Fig 2, only certain wing cells develop eyespot foci, while other wing cells do not develop any foci. The proposal of a mechanism that explains whether or not an eyespot focus forms in a given wing cell is one of the main objectives of the current study. We assume that the key determinant of focus point selection is in the overall venation system of the wing disc. Following Nijhout [9], we assume that veins act as sources of one of the two diffusing reactants. To investigate the selection mechanism, we assume a hypothetical venation system, where wing cells of the wing disc are rectangular (see Fig 3), although for completeness we also illustrate that our results are robust to perturbations from this rectangular geometry by considering wing cells with curved boundaries and varying width in Section 3.2. Under these assumptions, we investigate whether the nature of the proximal boundary condition can determine focus point selection, i.e., the number and locations of focus points in the entire wing disc.

Mathematical description of the model
For the sake of simplicity, we approximate the wing cells by rectangular domains and assume that all wing cells are of the same size (see Fig 3) and are independent [4]. Our results appear insensitive to small perturbations to this geometry. The only difference between wing cells is assumed to be in the source value of the activator at the proximal veins. We propose variable boundary conditions on the proximal boundary in the anterior-posterior direction of the wing disc. We assume that the entire wing disc comprises seven wing cells on which eyespot foci patterns can form. Following Nijhout [9], in the interior of each wing cell, we employ the activator-inhibitor reaction-diffusion model of Gierer-Meinhardt (G-M) [15] that describes focus point formation. Therefore, our model consists of several sets of coupled G-M models; each set is posed in a single rectangular wing cell. We remark that under this framework, the focus point formation occurs independently in each rectangular wing cell.
Let n seq denote the number of wing cells, typically seven. For the i-th (i = 1,. . .,n seq ) wing cell, the boundary conditions for the activator concentration (a 1 ) are Dirichlet (fixed) on the proximal boundary Γ p,i and the wing veins Γ v,i , Γ v,i+1 , and Neumann (zero flux) boundary conditions on the wing margin Γ m,i (i = 1,. . .,n seq ) (see Fig 3). The boundary conditions for the inhibitor concentration (a 2 ) are zero flux on all four boundaries of each rectangular wing cell. The Dirichlet boundary condition on each vein Γ v,i is the same for each vein. The initial conditions are taken to be the spatially homogeneous positive steady state solutions of the G-M equation. Thus our model for selection of focus points consists of n seq independent G-M equations. Let us denote by O i the i-th wing cell with boundaries, Γ m,i (wing margin), Γ v,i , Γ v,i+1 (veins) and Γ p,i (proximal boundary).
where the reaction function f ! ða ! Þ is given by with α, κ 1 , κ 2 , κ 3 > 0. This choice of reaction kinetics implies the existence of a positive steady state a ! ss of the ordinary differential equation (ODE) system and this is given by k 3 Other than the prescribed proximal boundary condition uðx ! Þ in Eq (3.1), each of the wing cells is assumed to be identical with the same source terms from the wing veins, diffusion coefficients and reactions. The boundary conditions at the veins are taken to be constant at twice the steady state of the activator, i.e.,ã ¼ 2 a 1 ss (following [9]).

Simulation results of the model with prescribed boundary conditions
We now present numerical simulations illustrating that the proximal boundary condition can act as a determinant of whether or not a focus point forms in a given wing cell. We use the finite element method, derived and analyzed in Lakkis et al. [16], for all the simulations approximating the equations on meshes with 33025 degrees of freedom and using a time-step of 10 −3 . We take the parameter values for the reaction kinetics and diffusion coefficients to be those given in Table 1. The majority of the results we report on remain qualitatively unchanged with small changes in the parameter values (10%). We consistently show only snapshots of the numerical solution of the activator concentration corresponding to Eq (3.1). The inhibitor concentration profile is in-phase and hence its snapshots are omitted. We start by considering prescribed boundary conditions on the proximal boundary, i.e., the function uðx ! Þ in Eq (3.1) is a given function. We consider the following three cases for the quadratic proximal boundary condition.

Constant boundary condition.
We first consider the case that the proximal boundary condition is constant in each wing cell, i.e., it is a piecewise constant discontinuous function over the whole proximal anterior-posterior boundary. In Fig 4, we show simulation results of Eq (3.1) on wing cells with constant proximal boundary condition of the form k p a 1 ss , where k p = 0, 1 and 2 (reading from left to right in each row) and a 1 ss is the (activator) steady state value. Each wing cell is taken to be a rectangular domain of length (proximal-wing margin) Table 1. Parameter values used for all the simulations of Eq (3.1). three and width (anterior-posterior) two. We observe the formation of activator peaks along the centerline of each wing cell (even those that do not eventually possess focus points) that is a characteristic of Nijhout's model [17] and is observed in experiments [1,2]. In wing cells with an activator concentration of less than 2a 1 ss on the proximal boundary, as the midline peak recedes it leaves behind focus points (two columns on the left in Fig 4), whilst for the cell with activator concentration on the proximal boundary equal to 2a 1 ss , (right hand rectangle) the , an extra focus point is formed which originates at the proximal boundary, migrates to the interior of the wing cell and persists at the steady state. Thus, for this choice of parameter values and domain geometry, the number of focus points at steady state does not depend monotonically on the proximal boundary condition. As each wing cell only differs in terms of the proximal boundary condition, we see that the changes in the proximal boundary condition can act as a determinant of focus point formation. The piecewise constant boundary profiles considered so far are only an approximation and it is likely that the real activator boundary profile may appear as a continuous smooth function.
3.2.2 Concave and convex boundary conditions. We consider the following two additional proximal boundary condition profiles, a concave profile: and a convex profile: where w is the width of the wing cell. Fig 5 shows simulation results of Eq (3.1) on wing cells together with two profiles of the proximal boundary conditions given by Eqs (3.3) and (3.4). As previously, each wing cell is taken to be a rectangular domain of length (proximal-wing margin) three and width (anteriorposterior) two. We once again observe the formation of activator peaks along the centerline of each wing cell and as this midline peak recedes, it leaves behind a spot in the wing cell with the concave boundary condition whilst with the convex boundary condition the peak completely recedes leaving behind no spot. We have performed a number of other simulations (results not shown) with spatially varying (within each wing cell) proximal boundary conditions and we observe analogous behavior to this simulation, namely that the value of the boundary condition in the middle of the proximal boundary of the wing cell is a key in determining whether or not a focus point forms. We further note that by appropriately super imposing together boundary conditions as in As a robustness test of sensitivity of the observed behavior to the geometry, we relax the assumption of rectangular geometries and work on a geometry closer to that of the real wing cells shown in Fig 2. We consider a wing cell whose width increases as we move in the proximal-distal direction and we consider curved proximal and wing margin boundaries. The specific geometry for which we present results is defined by the following boundaries: Γ v,1 is taken to be the line between the points (-0.8, 3) and (-1, 0), Γ v,2 is taken to be the line between the points (0.83, 3) and (1, 0) and the proximal and wing margin boundaries are taken to be curves given by shows simulation results of Eq (3.1) on wing cells with curved proximal and wing margin boundaries and increasing width towards the wing margin. We observe analogous behaviors to the rectangular domain case with larger values of the Dirichlet proximal boundary condition inhibiting the formation of a focus point. In the case of values at the steady state for the proximal Dirichlet boundary condition, we note that in contrast to the rectangular domain case only one focus point is observed (results not included in the interests of space), with the steady state profile similar to those of the zero Dirichlet boundary condition (Fig 6 (left)).
Finally, we investigate the dependence of the number of focus points on the aspect ratio of the wing cell (see figures (a)-(h) in S1 Appendix for detailed numerical simulations). Figures  (g) and (h) in S1 Appendix show steady states of the simulation of Eq (3.1) on wing cells with constant proximal boundary condition at zero and 2 times the steady state value. The length (proximal-distal) is held fixed at 3 and the width (anterior-posterior) is varied between 1.5 and 3 (i.e., the aspect ratio varies from 2 to 1). We observe a monotonic dependence of the number of focus points on the aspect ratio. We have also simulated cases (results not shown) where the proximal boundary condition is asymmetric and where the boundary condition along the veins is varied, rather than the fixed Dirichlet conditions presented above. Although it is possible through careful tuning of the parameter values and boundary conditions to generate focus points which are not circular (spots) such as arc shaped foci as well as focus points which are positioned away from the centerline of the wing cell, the predominantly observed behavior is the generation of circular spot shaped foci positioned along the centerline of the wing cell at steady state.

Simulations of variations in focus point patterning observed in nature
We now present some simulation results together with experimental images of real specimens, which illustrate the capability of the model to describe naturally occurring variations in eyespot focus point patterning.
3.3.1 Development of focus points in the wing disc during eyespot determination.  the interior of the rectangle. For this simulation we selected boundary conditions to be four times the steady state on the completely developed veins and the proximal boundary and twice the steady state on the incompletely developed vein. The interior boundary was modeled as a Dirichlet boundary only for the activator whilst for the inhibitor all the boundaries were taken to be zero-flux. The inclusion of such an interior boundary in the finite element simulations is straightforward once a triangulation is defined over the desired geometry.
We see that with this choice of boundary profiles the resulting focus point distribution is similar to that observed in the case of abnormal wing venation. The results suggest that the incomplete vein may constitute a smaller source of activator than completely developed veins and this could account for the variation in the position of focus points and hence the resultant eyespot pattern.
3.3.3 One eyespot splits into two eyespots through the addition of a vein. Finally, we conclude this subsection with another example of an abnormal eyespot pattern on the hind wing of the butterfly Y. arugus, which at first glance appears incompatible with the current incompatible with our model, as the results in S1 Appendix (i.e., the influence of aspect ratio on focus point determination) show that, in general, reducing the width of the wing cell leads to the formation of fewer foci. However once again if we assume that the abnormal case corresponds to a change in the venation system, specifically, a change in the boundary conditions at the newly formed vein then the model is capable of generating results consistent with the experimental observations. Fig 9(C) shows results of a simulation on two rectangles of length three and width one, i.e., wing cells of half the usual width. The black line in the figure indicates the new vein. We assume that this new vein acts as a homogeneous Dirichlet boundary for the activator. The proximal boundary is also taken to be a homogeneous Dirichlet boundary for the activator whilst the other pre-existing vein is assumed to be a Dirichlet boundary, with the concentration at the boundary at twice the steady state. Due to symmetry, we only solve on a single wing cell and simply reflect along the line of the new vein. We clearly see the formation of two focus points, one in each thin wing cell. Hence under our model one prediction is that for abnormal scenarios such as those in Fig 9 lower activator source strengths on the abnormally inserted veins may account for the abnormal patterning.

A 2-Stage Model for Focus Point Selection
In Section 3.2, we illustrated that the consideration of different proximal boundary conditions is sufficient to explain focus point selection in a single wing cell. In order to present a complete model for focus point selection, it remains to develop a mechanism for the generation of the proximal boundary profiles. We propose a 2-stage process whereby the first stage consists of the formation of the pattern generating the proximal boundary profiles and the second stage consists of the focus point formation model described in Section 3.1. Although we consider a 2-stage model in this work for simplicity, however, it is certainly of interest mathematically and may be biologically important to consider models where the boundary profile pattern formation process occurs on the same timescale as the focus point formation process. Such a coupled bulk-surface system may be an interesting direction for future research.

A model for the generation of the proximal boundary condition and a 2-stage model for focus point patterning
For the generation of the proximal boundary profiles, we propose an 1-dimensional (1D) pattern formation model posed on the proximal boundary Γ p = S i Γ p,i , the union of the proximal boundaries of the wing cells. Clearly, a large variety of models could generate boundary conditions of the form considered in the previous section, the 1D model we present here is just one concrete example.
To illustrate the modeling, we work with a concrete example, the activator depleted substrate model of Schnakenberg [19] (see also Murray [20]). The reaction-diffusion system (RDS) is posed on the anterior-posterior margin of the entire wing disc, i.e., the proximal boundary Γ p and we assume zero-flux boundary conditions. We consider the following dimensionless RDS for the concentrations of two chemicals (activator and substrate): Find u ! ðx; tÞ ¼ ðu 1 ðx; tÞ; u 2 ðx; tÞÞ T such that where, d 1 , d 2 , k 1 and k 2 are all positive constants. r Γ and Δ Γ denote the surface gradient and Laplace-Beltrami operators, respectively. Usually the function γ appearing in Eq (4.1) is taken to be a positive constant (which may be interpreted as being related to the domain size or alternatively may be interpreted as a reaction rate [20]). However, in general due to the inherent uniform wavelength associated with Turing patterns, we believe it is not possible to generate boundary profiles such as those considered in Section 3.1, which allow focus points to be generated in arbitrary wing cells with constant parameters if one considers only two component RDSs. We propose here the 1 D continuous two component RDSs (4.1) with a spatially varying γ to obtain different proximal boundary profiles. We note that 3 or more component RDSs have much richer behavior than 2 component RDSs [21 and 22] and hence may be attractive candidates for the generation of the anterior-posterior pattern.We remark, that this anteriorposterior patterning may occur on a different timescale to the focus point formation process, and it may even occur at an earlier stage of the focus point development, and it would then lay down a pre-pattern for the formation of the proximal boundary profile. We now describe the 2-stage model we propose for the modelling of focus point patterning on butterfly wings: Stage 1: In the first stage, the 1D RDS (4.1) is solved on the proximal boundary Γ p (i.e., the union of the proximal boundaries of the wing cells) to steady state. In the next section, we present firstly (1) simulation results of Eq (4.1) in cases where the function γ is constant, and secondly (2) simulation results of Eq (4.1) with a spatially varying γ. We reused the parameter values given in Table 1 for the system (3.1) and the remaining parameters were taken as shown below. Thus as the parameter c p,2 = 0, the Dirichlet boundary conditions u in Eq (3.1) are simply given by one third of the (patterned) steady state activator concentration, u 1 of the 1D RDS (4.1). The initial conditions for the 1D RDS are taken as small quasi-random perturbations around the uniform steady state ðk 1 þ k 2 ; k 2 = ðk 1 þ k 2 Þ 2 Þ T and the initial conditions for the bulk RDS are taken as in Section 3.1 (uniform steady state values). In all the simulations, we assume the idealized geometry depicted in In this case, γ is below the critical value for the onset of diffusion driven-instability globally and the solution to the 1D RDS (4.1) simply converges to the uniform steady state. Hence, we have a constant value for the proximal boundary condition. As the choice of the coupling coefficient c p,1 is such that the proximal boundary condition ðc p;1 u 1 Þ is below the critical value for focus point formation, we generate a focus point in every wing cell. The resultant pattern is similar to that observed in the developing wing disc of B. anynana as shown in Fig 10 (right hand column).

Simulation results of the 2-stage model
Case 2: No focus points on the wing: Fig 11 shows results of the coupled model with γ(x) = 1. In this case, γ is above the critical value for the onset of diffusion driven instability and we obtained a patterned steady state solution to the 1D RDS (4.1). The solution profile is one of seven equally spaced activator peaks in the domain with the activator profile on each proximal boundary ([0,2], [2,4],. . .) appearing similar in shape to the convex profile of Fig 5. When we simulate the 2-stage model with this boundary profile and the selected coupling coefficient c p,1 , we observe similar behavior to the simulations of Fig 5 (with the convex profile) with no focus points forming in any of the wing cell.
Beldade et al. [4] did artificial selection to examine how the relative size of the anterior and posterior eyespots on the dorsal forewing of B. anynana can be changed in a laboratory Case 3: Two focus points on the dorsal hindwing of Precis coenia: Brakefield et al. [2] showed that the fifth-instar Precis coenia hindwing imaginal disc exhibits two spots of Dll expression, which correspond to the future positions of two eyespots on the adult hindwing. In Fig 12, the top right photo (b) shows the adult P.coenia dorsal hindwing with two eyespot focus points, and the bottom right (d) is the fifth instar P. coenia hindwing imaginal disc displaying a pre-pattern with two focus points of Dll expression, which correspond to the position of the eyespots on the adult dorsal hindwing.
By selecting, the parameters in the two component RDS posed on the proximal boundary corresponding to the third case described above, we generate patterns with a larger wavelength (due to the increased diffusion coefficients). As shown in Fig 12(A), the resultant steady state consists of an activator pattern with only two interior minima in the domain corresponding to the proximal boundaries of wing cells two and five. When we simulate the 2-stage model with this boundary profile and the selected coupling coefficient c p,1 , we observe the formation of foci in wing cells 2 and 5 and no foci in the other wing cells similar to the experimental observations of P. coenia as shown in Fig 12 (right hand column).
4.2.2 Simulation results of (4.1) with a spatially varying γ in the anterior-posterior direction. As mentioned in Section 4.1, through the consideration of 2-component RDSs with constant parameters, for the 1D patterning mechanism appears insufficient to generate boundary profiles leading to focus points in an arbitrary wing cells. A major difficulty lies in the fact that Turing patterns typically possess a constant wavelength over the domain, hence patterned profiles in one region of the domain with no patterning (convergence to the homogeneous steady state) in another region is not possible. However, such a pattern distribution is easily achieved through the consideration of systems with spatially varying parameters. See for example [23] for previous work in this direction. In upper figures in S2 Appendix, we present simulation results of Eq (4.1) (as part of the coupled model) that illustrate that a system of the form (4.1) with a spatially varying γ, specifically an anterior-posterior gradient in γ can restrict patterning to certain portions of the wing (see also the corresponding bottom figures in S2 Appendix). RDSs with spatially varying parameters have been the subject of much study in the literature, for example [24]. We remark that numerical studies suggest patterning can be restricted to certain portions of the domain through the use of spatially varying parameters similar to the results we report on the current work. In [25,26], the authors model the regulation of digit patterning of developing vertebrate limb buds by Hox genes using a Turing RDS. Their results indicate that changing the kinetic parameters can influence the wavelength of the resultant pattern. They also consider spatial gradients in kinetic parameters, and show that this allows the robust formation of a striped pattern with a given orientation that models digit formation.

Summary and Discussion
In this study, we presented a model for the selection and distribution of eyespot focus points on the wings of Lepidoptera. The basic idea of the model is that the wing cells, in which eyespot foci are formed, are selected by the source value of an activator on the proximal veins of the entire wing disc. Specifically a variable proximal boundary condition in the anterior-posterior direction of the entire wing disc governs focus point selection. Through numerical simulations on idealized wing disc geometries, we illustrated that this proximal wave-like boundary condition can govern the number and locations of eyespot focus points on the wing surface. As a result, the model could provide a plausible mechanism for the selection of global eyespot focus points on wing discs of some butterfly species such as Junonia (or Precis) coenia, Bicyclus anynana, and Ypthima arugus in the Nympalidae family.
Our study suggests that a key factor that determines focus point selection could be in the overall venation system of the wing disc. We assumed that the veins are sources of the activator [9], i.e., they act as Dirichlet boundaries in the mathematical model and numerical simulations. We first considered a number of different prescribed boundary conditions on the proximal boundary, and our results show that one may construct boundary profiles (which could be smooth and continuous or discontinuous across the whole wing disc) such that focus points may be selected in any wing cell and hence may reproduce the variety of focus point distributions observed in experiments. To complete the model, we then proposed a simple 1D reaction-diffusion model for the generation of the proximal boundary condition profile, that is, a surface Turing system posed in the anterior-posterior direction on the proximal vein.
We stress that the key factor is a change in source values of the activator on the proximal veins in the anterior-posterior direction of the wing disc and this change may be realized by a variety of different patterning mechanisms. One of our main results is that under our model, wing cells in which eyespot focus points are generated need to have lower source strength of the activator on the proximal boundary than wing cells that do not produce focus points. We also stress that the number and locations, that is, the global distribution of eyespot foci on the wing disc could be largely controlled by two gradients along two different directions, that is, the first one is the gradient in spatially varying parameters such as the reaction rate γ in Eq (4.1) along the anterior-posterior direction on the proximal boundary of the wing cells (see Section 4.2.2 and S2 Appendix), and the second one is the gradient in source values of the activator along the veins in the proximal-distal direction of the wing cell. The first gradient could determine the number and locations of foci on the wing cells. The second gradient could control the location of the focus point in the proximal-distal direction within the wing cell (see also discussion in Section 3.2).
The role of boundary conditions in the determination of patterning generated by reactiondiffusion systems has been the subject of previous work. For example, Dillon et al. [27] show that changes in boundary conditions can have a profound influence on the solutions both in terms of existence and uniqueness and in terms of the stability of patterns. Page et al. [24] investigate Turing systems with a discontinuity in a kinetic parameter and show that such a system may be decomposed into systems with constant parameters and anomalous boundary conditions. They further show that such systems may exhibit spatial patterns outside the classical Turing space. Barrio et al., [28] and Aragon et al., [29] consider the role of boundary conditions in Turing reaction-diffusion system models for pigment patterns on the skin of fish. In particular, they observe that the consideration of different boundary conditions may increase the number of scenarios such models are capable of explaining. In light of our numerical results and the theoretical and numerical works mentioned above, further mathematical investigations in the same direction as the works above into the role of boundary conditions in patterning by Turing systems are certainly warranted.
The current model framework consists of wing discs that are the union of several (identical) rectangular wing cells, although we also considered a somewhat more representative wing cell with curved (proximal and wing margin) boundaries and varying width. This modeling framework might be improved by considering a more realistic geometry of the wing cells, rather than the rectangular cells considered in this study, although experimental results suggest that during the time at which focus points form, a rectangular cell is a good approximation. Preliminary numerical results suggest that other forms of boundary conditions can generate more complicated foci such as arc shaped foci or multiple foci of different sizes in a wing cell. Such models could perhaps reproduce some of the diversity of focal shapes that are occasionally observed in nature. The exploration of such possibilities are reserved for future work. Here our objective has been to provide a proof of concept that anterior-posterior patterning alone may determine focus point selection.
Supporting Information S1 Appendix. The influence of aspect ratio on eyespot focus point determination. The figures in S1 Appendix show snapshots of the activator concentration corresponding to the solution of Eq (3.1). The wing cells are taken to be rectangular and of fixed length equal to three, but the width is now varied with the width taken to be 1, 1.5, 2, 2.5 and 3 reading from left to right in each subfigure. In the left hand column, the boundary conditions for the activator on the proximal boundary (top) of each rectangular cell are taken to be zero, and in the right hand column, they are set to twice the steady state value. Initially in all but the thinnest wing cell (left hand most in each column), a vertical stripe of high activator concentration is generated originating from the zero-flux distal boundary (bottom). As the width of the wing cells is increased, the midline peak starts to generate multiple spots as it recedes with insertion of new spots in the proximal-wing margin, and in the anterior-posterior direction, both exhibited. There appears to be a monotonic relationship between aspect ratio and number of focus points with wider wing cells (with a fixed length) exhibiting more focus points at steady state. (TIF) S2 Appendix. Restricting focus point formation to regions of the wing through the use of spatially varying parameters in the anterior-posterior 1D RDS. As mentioned in Section 4.2.2, 2-component RDSs with constant parameters alone for the 1D patterning mechanism appear insufficient to generate boundary profiles leading to focus points in arbitrary wing cells. However, such a pattern distribution is achieved through the consideration of systems with spatially varying parameters such as the reaction rate γ(x) in Eq (4.1). To illustrate this effect, in figures in S2 Appendix, we report on the steady states for the 1D systems obtained using a monotonically increasing gradient for the reaction rate γ(x) = (x/14) p , with p = 1, 2 and 8 ((a), (b), and (c), respectively). The remaining parameters were taken to be with κ 1 = 0.1, κ 2 = 0.9, d 1 = 0.01, d 2 = 1, c p,1 = 1/3 and c p,2 = 0. As the gradient of this function is decreased (smaller p), focus points form only closer to the anterior margin whilst for larger gradients focus points can be made to form on almost the entire wing. The case p = 0, corresponds to Fig 11 with no focus points forming, whilst formally in the limit p!1, a focus point forms in each wing cell, as in Fig 10, as the 1D RDS solution would be close to the steady state value due to the initial conditions. (TIF)