## Mathematical modeling of the eyespots in butterfly wings

1

Kang-Ling Liao^{1,2,∗} Wei-Chen Chang^{3}, Jeffrey M. Marcus^{2}, Jenn-Nan Wang^{4}

2

1 Department of Mathematics, University of Manitoba, Manitoba, R3T 2N2, Canada

3

2 Department of Biological Sciences, University of Manitoba, Manitoba, R3T 2N2, Canada

4

3 Department of Mathematics, National Taiwan University, Taiwan, R.O.C.

5

4 Institute of Applied Mathematical Sciences, National Taiwan University, Taiwan, R.O.C.

6

∗ Corresponding author. E-mail address: Kang-Ling.Liao@umanitoba.ca

7

## Abstract

8

Butterfly wing color patterns are a representative model system for studying biological pattern formation,

9

due to their two-dimensional simple structural and high inter- and intra-specific variabilities. Moreover,

10

butterfly color patterns have demonstrated roles in mate choice, thermoregulation, and predator avoidance

11

via disruptive coloration, attack deflection, aposematism, mimicry, and masquerade. Because of the

12

importance of color patterns to many aspects of butterfly biology and their apparent tractability for

13

study, color patterns have been the subjects of many attempts to model their development. Early

14

attempts focused on generalized mechanisms of pattern formation such as reaction-diffusion, diffusion

15

gradient, lateral inhibition, and threshold responses, without reference to any specific gene products. As

16

candidate genes with expression patterns that resembled incipient color patterns were identified, genetic

17

regulatory networks were proposed for color pattern formation based on gene functions inferred fromother

18

insects with wings, such as Drosophila. Particularly detailed networks incorporating the gene products,

19

Distal-less (Dll), Engrailed (En), Hedgehog (Hh), Cubitus interruptus (Ci), Transforming growth factor-β

20

(TGF-β), and Wingless (Wg), have been proposed for butterfly border ocelli (eyespots) which helps the

21

investigation of the formation of these patterns. Thus, in this work, we develop a mathematical model

22

including the gene products En, Hh, Ci, TGF-β, and Wg to mimic and investigate the eyespot formation

23

in butterflies. Our simulations show that the level of En has peaks in the inner and outer rings and

24

the level of Ci has peaks in the inner and middle rings. The interactions among these peaks activate

25

precursor cells of pigments to generate white, black, and yellow pigments in the inner, middle, and

26

outer rings, respectively, which captures the eyespot pattern of wild type(Bicyclus anynana)butterflies.

27

Additionally, our simulations suggest that lack of En generates a single black spot and lack of Hh or Ci

28

generates a single white spot, and adeficiency of TGF-β or Wg will cause the loss of the outer yellow

29

ring. These deficient patterns are similar to those observed in the eyespots of Vanessa atalanta, Vanessa

30

altissima, and Chlosyne nycteis. Thus, our model also provides a hypothesis to explain the mechanism

31

of generating the deficient patterns in these species.

32

Keyword: Butterfly wings, pattern formation, reaction-diffusion model, gene expression.

33

## Introduction

34

Butterfly wing color patterns are an attractive model system for studying biological pattern formation.

35

Color patterns are particularly suitable for such studies because they are structurally simple and two-

36

dimensional, they consist of clearly defined subunits, and they are highly inter- and intra-specifically

37

variable [1, 6, 14, 33, 39, 47]. Butterfly color patterns have demonstrated roles in mate choice [39, 59],

38

thermoregulation [24], and predator avoidance (including camouflage [68], disruptive coloration [57],

39

attack deflection [21], aposematism [15], mimicry [58], and masquerade [65]). Because of the importance

40

of color patterns to many aspects of butterfly biology and their apparent tractability for study, color

41

patterns have been the subjects of many attempts to model their development.

42

Each wing surface consists of a flat and static monolayer of epidermal cells [50]. A subset of the

43

epidermal cells differentiate into scale cells [62, 64] which will synthesize pigments [26, 28] to generate

44

color patterns including border ocelli (also known as eyesopts) [50]. The number, location, and size of

45

the eyespots differ among species of butterflies [50]. The position and shape of the wing color pattern

46

are determined by the locations of signalling sources from the wing veins and wing margin [28, 47, 50].

47

Additionally, the pattern formation requires a two-step process: the determination of the distribution

48

of discrete signalling sources for the color pattern during the last larval instar and the differentiation of

49

the surrounding pattern during the pupal state [46]. The wing patterns (i.e., the background pattern

50

or global pattern) are composed by five pattern elements: i) ripple patterns are the rhythmical patterns

51

covering the whole wing surface, ii) dependent patterns are the pattern depended on the lacunae of the

52

pupal wing, iii) crossbands are the band pattern alone the anterior to the posterior margin of the wing,

53

iv) eyespots are the pattern consists of concentric rings with different colors, and the v) color fields are

54

the large areas of the wing surface with color [44]. The position, number, size and color of eyespots

55

are determined in a developmental pathway that is independent of other pattern elements and body

56

structures [8]. The eyespot appears from an inductive organizing center, the focus,which isa signalling

57

source of a morphogen to determine the pigment of surrounding cells [8]. The color of the wing patterns

58

is determined by the pigment generated by scale cell surface features that reflect light [46, 61, 62, 64].

59

In [8], Brakefield et al. defined the developmental pathway for eyespot formation and then investigated

60

how these pathways affect the numbers and sizes of eyespots in the squinting bush brown butterfly,

61

Bicyclus anynana (Lepidoptera: Nymphalidae: Satyrinae). The experiments inBicyclus anynanashowed

62

that the dynamics of the expression of Distal-less (Dll) gene can be used to categorize the eyespot

63

formation into the following four stages [49]. Stage I - The larval prepattern: In the larva, a high level

64

expression of Dll protein appears as a broad band and stripes down the middle of each wing subdivision to

65

creates the potential focal pattern. Stage II - The focal determination: The Dll protein accumulates and

66

is stabilized at the tips of these stripes and then diffuses to form stable circular spots of Dll expression.

67

Stage III - The focal signaling: In the pupa, the high level of Dll expression expands to a broader circular

68

region where the region is determined by the signaling from the epithelial focus. Stage IV - Differentiation:

69

The positions of these spots of Dll expression become the central regions of the eyespots. Some graded

70

morphogen appears across the radius of the eyespot region such that the surrounding cells of foci generate

71

different pigmented scale types, according to the level or type of signal they receive and their location

72

within the wing [8]. Based on this eyespot formation process, the Dll expression can be used to determine

73

the position and number of foci and detect signaling from the focus [8]. Within the wing epidermis,

74

the signalling molecules move through the extracellular medium by diffused through gap junctions [46].

75

Additionally, there is no cell migration within the wing epithelium, so the pattern is mainly affected by

76

the cell differentiation in responsive to chemical signalling, instead of the cell movement [46, 50].

77

Many studies employ the activator-inhibitor type of models used by Turing [66] to investigate the

78

wing pattern formation of butterfly [5, 13, 32, 43–45, 47, 48, 50, 61]. The first proposed model type is the

79

gradient model that the morphogens produced by central cells of eyespots diffuse to the surrounding

80

cells, and then the surrounding cells differentiate into discrete rings based on the received morphogen

81

concentration [39, 44, 45, 47]. Nijhout provided gradient models based on the distribution information

82

from focus of ocellus [44] to show that the foci are the sources of a diffusing chemical to activate color-

83

specific biosynthetic pathways [45, 47] and to show that the chemical signal depending on cell position

84

generates the surrounding patterns [47]. Nijhout also proposed a model involving the signaling from an

85

activator in a lateral inhibition reaction to generate the required spatial distribution of sources and sinks

86

such that the model can generate the eyespot patterns [46, 48]. Murray also proposed a diffusion model

87

incorporating a diffusing morphogen resulting in the activation of a biochemical gene to generate the

88

wing patterns in lepidopteran [43].

89

Bard et al. [5] provided a diffusion equation incorporating morphogen sources at the foci and sinks

90

at the wing margin with appropriate diffusion throughout the wing. Since the morphogen concentration

91

determines the pigments generated by scale cells, this model is able to generate wing patterns for different

92

species of butterflies [5]. In [61], Sekimura et al. created a modified Turing mechanism reaction-diffusion

93

model, involving different regions in the wing and spatially dependent morphogen, to generate the global

94

pattern on the wing of P. dardanus on a geometrically accurate wing domain. Sekimura et al. used this

95

model to investigate the parameter values for mode selection, threshold values for color determination,

96

wing shape and boundary conditions [61], and then predict the global effect on wing patterns in cutting

97

experiments [32] which cannot be generated from the model mentioned in [47]. In [13], Dilao and Sainhas

98

provided a reaction-diffusion model, involving two diffusive morphogens for the first eyespot ring for-

99

mation and the modification of wing background pigment precursors, to generate the general structural

100

organization of eyespots.

101

These activator-inhibitor models provide analysis of color pattern formation based on the interaction

102

between generalized activator and inhibitor morphogens, but they still lack detailed information concern-

103

ing the specific morphogens and signaling processes involved in the development of butterfly eyespots.

104

Thus, based on the findings from activator-inhibitor models, researchers started to build models with

105

detailed structural analysis of eyespots. In [52–54], Otaki provided a simple uniformly decelerated mo-

106

tion model describing the interaction between the morphogenic signals and parafocal elements (PFEs) to

107

investigate the universally morphological feature inside-wide rule of eyespots: one eyespot contains one

108

inner core black ring and an outer black ring. In [63], Sekimura et al. built a spatially two-dimensional

109

reaction-diffusion system model with non-homogeneous Dirichlet boundary conditions to investigate the

110

mechanism for determination of the number and location of eyespots. Their simulation results suggested

111

that the morphogen concentration along the proximal vein is the main factor to control the distribution

112

of eyespots and this observation is robust to the proximal boundary condition [63].

113

In [14], Evans and Marcus provided several reaction-diffusion models to generate the concentrations of

114

gene expression during the eyespot formation inBicyclus anynana and Junonia coenia. Comparing their

115

simulation result with experimental data for these gene expressions, they made the following conjecture.

116

In the eyespot foci (i.e., the inner ring of the eyespot), firstly, the expression of Notch gene induces the

117

co-expression of Notch and Distal-less (Dll). Secondly, the expression of Dll activates the expression of

118

the gene Engrailed (En). Thirdly, the expression of En shows positive associations with the expression

119

patterns of hedgehog (Hh) transcript. Fourthly, the expression of hh triggers the expression of Patched

120

(Ptc) and the transcription factor Cubitus interruptus (Ci), which suggested that the expression of

121

hh promotes the expression of Ci. However, the hh expression inhibits the production of Ci within

122

intracellular reaction suggesting that the intracellular function of hh actually suppresses the expression

123

of Ci. In [34], Marcus and Evans used the model mentioned in [14] to generate the wing patterns in

124

two mutants, the comet mutation with a series of comet-shaped eyespot foci and the Cyclops mutation

125

with failure in wing vein formation. In [33], Marcus provided more complete information about the

126

downstream pathways of En gene expression, such as how the interaction between the TGF-β signaling

127

and Wingless (Wg) signaling affects the localization of eyespots. Additionally, Marcus also explained the

128

mechanism for generating different pigments in different region within the eyespots: the concentration of

129

En in the inner and outer rings is used to generate the white and yellow pigments in the inner and outer

130

rings respectively, and the concentration of Ci is used to active the generation of black pigment in the

131

middle region of the eyespots [33].

132

A summary of these mathematical models and the focus of this workarelisted in the Table 1. In Table

133

1, all model types (i.e., the second column)incorporate cellular responses to diffusing morphogens. Mod-

134

els that incorporated generalized reaction-diffusion interactions are labeled as reaction-diffusion models.

135

Gradient models incorporate a diffusion term that represents the gradient of a single activator, while Tur-

136

ing reaction-diffusion models include two morphogens: one activator and one inhibitor. Therefore, these

137

models used different mathematical approaches and cannot be mutually replaced. We use the Wolpert

138

Positional Information (PI) theory mentioned in [69, 70] to define different types of patterning mecha-

139

nisms: In the fourth column, the de-novo patterning refers to the initial specification of developmental

140

pre-patterns and organizers whereas the fine-scale patterning represents the effect of those pre-patterns

141

and organizers on the surrounding cells and tissues.

142

Table 1. Summary of mathematical models for eyespot formation. The reference, type, and the stage defined in [8] of the models are listed in the first, second, and third columns, respectively. The patterning type defined in [69, 70] is shown in the fourth column. The fifth column shows whether the model incorporates molecular or genetic information. The last column describes whether the simulation results are only hypothetical or have support from experiments.

Model Type Stage De-novo v.s. fine-scale Molecular & Hypothesis v.s.

patterning genetic details experiment

Nijhout 1978 Gradient model II De-novo Excluded Experiment

[44] [16, 44]

Murray 1981 Turing reaction-diffusion model II & III De-novo/Fine-scale Excluded Experiment

[43] (inhibitor-activator) [43]

Bard & French 1984 Turing reaction-diffusion model II & III De-novo/Fine-scale Excluded Hypothesis

[5] (inhibitor-activator)

Nijhout 1980, Nijhout 1991 Gradient model II & III De-novo Excluded Experiment

[45, 47] [17, 45]

Nijhout 1990, Nijhout 1994 Turing reaction-diffusion model I &II De-novo Excluded Experiment

[46, 48] (inhibitor-activator) [63]

Sekimura et al 2000 Turing reaction-diffusion model I De-novo Excluded Experiment

[61] (inhibitor-activator) [32]

Dila˜o& Sainhas 2004 Turing reaction-diffusion model IV Fine-scale Included Experiment

[13] (inhibitor-activator) [17, 48]

Evans & Marcus 2006 Reaction-diffusion model II De-novo Included Experiment

[14] [14, 23, 27]

Marcus & Evans 2008 Reaction-diffusion model II De-novo Included Experiment

[34] [8]

Otaki 2011, Otaki 2012 Gradient model IV Fine-scale Excluded Experiment

[52–54] [54]

Sekimura et al 2015 Turing reaction-diffusion model I & II De-novo Excluded Experiment

[63] (inhibitor-activator) [63]

This work Reaction-diffusion model III Fine-scale Included Hypothesis

In the current work, based on the interaction amount En, Hh, Ci, TGF-β, and Wg gene productsin

143

Bicyclus anynana and Junonia coenia mentioned in [14,33,34], we develop a system of partial differential

144

equations (PDEs), including the concentrations of En, Hh, Ci, TGF-β, and Wg proteins, to mimic the

145

eyespot formation in Bicyclus anynana butterflies. Our simulation shows that the concentration of En

146

has one peak in the inner ring and one peak in the outer ring, as well as the concentration of Ci has one

147

peak in the inner ring and one peak in the middle ring. High concentrations of En and Ci in the inner

148

ring triggercells to express the biosynthetic pathway responsible for making white pigment (pteridine),

149

high concentration of Ci in the middle ring activatesce;;s tp express the biosynthetic pathway responsible

150

for making black pigment (dopa melanin), and the high concentration of En in the outer ring initiates

151

expression of components of the biosynthetic pathway for the yellow pigment (pheomelanin). Therefore,

152

our simulation captures the generation dynamics of white, black, and yellow pigments in the inner, middle,

153

and outer rings respectively which fits the pattern of wild type Bicyclus anynana butterflies described

154

in [33]. On the other hand, we use this model to predict the eyespot patterns in knockout mutants. Our

155

simulations display three types of degenerated patterns: (i) a single black spot observed from Vanessa

156

atalanta is caused by the deficiency of En; (ii) a single white spot observed from Vanessa atalanta and

157

Vanessa altissima is caused by the deficiency of Hh or Ci; and (iii) loss of the outer yellow ring shown

158

from Chlosyne nysteis is caused by the loss of TGF-β or Wg. A summary of our work for the eyespot

159

patterns in wild type Bicyclus anynana butterflies and these null mutants is shown in Fig 1. Finally,

160

our sensitivity analysis shows that (i) increasing the production rate of Ci in inner ring, or reducing the

161

production rate of Hh in inner ring or the production rate of Ci in middle ring promotes the white pigment

162

formation in the inner ring; (ii) enhancing the production rate of Ci in middle ring or degradation rate

163

of En, or reducing the production rates of En, Hh, or Ci in inner ring, or the diffusion rate of TGF-β

164

promotes the black pigment formation in the middle ring; and (iii) increasing the production rate of En

165

in outer or inner rings promotes the yellow pigment formation in the outer ring.

166

As listed in Table 1, most pervious mathematical models are based on theTuringreaction-diffusion

167

model, and hence lack the molecular and genetic information and mainly capture the behavior in stage II

168

during the eyespot pattern formationwhen the eyespot focus is specified. Thus, these models cannot be

169

used to investigate how the gene expression results in eyespot ring formation in stages III and IV. Here

170

we attempt to incorporate these missing stagesin our model, so the main contribution of our work is to

171

provide a different modeling approach incorporating molecular information to investigate how the cells

172

react related tomorphogens according totheir distancefrom the focus in stage III.

173

## Results

174

### Mathematical model

175

Our mathematical model is based on the network described in Fig. 2. The variables that will be used

176

are listed below. The values of parameters are listed in Table 2 and are estimated by using experimental

177

data in the Method section:

178

E(x, t) = concentration of En protein at location x and time t with unit kD/cm, H(x, t) = concentration of Hh protein at location x and time t with unit kD/cm, C(x, t) = concentration of Ci protein at location x and time t with unit kD/cm, T (x, t) = concentration of TGF-β protein at location x and time t with unit kD/cm, W (x, t) = concentration of Wg protein at location x and time t with unit kD/cm.

An eyespot consists of inner, middle, and outer rings. Thus, it is reasonable to consider the radially

179

symmetric solutions and reduce the two dimensional spatial variable x to the one dimensional variable,

180

r, representing the radius of the eyespot. The total eyespot area in wild typeBicyclus anynana is around

181

11.2 mm^{2} [55] and the radius of the eyespot is around R0= 0.094 cm. For simplicity, we focus on the

182

eyespot region and nondimensionalize the distance between the center and boundary of the eyespot to

183

be 1 by rescaling R0(= 0.094 cm). In this work, we are mainly interested in the qualitative description

184

of the eyespot patterns. Hence, we explicitly set the inner, middle, and outer rings as follows

185

Ω_{in}:= {0 ≤ r ≤ 0.3}, Ω_{mid}:= {0.3 ≤ r ≤ 0.6}, Ω_{out} := {0.6 ≤ r ≤ 1}. (1)
Next, we define the functions χ_{in}(r), χ_{mid}(r), and χ_{out}(r) as

186

χ_{in}(r) = 1

1 + e^{30(r−0.15)}, (2)

χmid(r) = 1.5 × 1

1 + e^{−30(r−0.4)}× 1

1 + e^{30(r−0.5)}, (3)

χ_{out}(r) = 1

1 + e^{−30(r−0.8)}, (4)

to restrict reactions in the inner, middle and outer rings, respectively. Notice that the simulation results for wild type and null mutants butterflies will not be changed, if the regions Ωin, Ωmid, and Ωout are replaced by

Ωin:= {0 ≤ r ≤ a}, Ωmid:= {a ≤ r ≤ b}, Ωout:= {b ≤ r ≤ 1},

with 0 < a < b < 1 and the definitions of χin(r), χmid(r), and χout(r) are adjusted accordingly.

187

Figure 1. Summary of the relation between the system network and predicted eyespot pattern. (A) shows the system network among the five keygene expressions: En, Hh, Ci, TGF-β, and Wg of wild type Bicyclus anynana,based on experimental evidences. The detailed explanation of the network is shown in Fig. 2. (B) provides the gene expression profiles at 16 hours over the radius of the eyespot in wild type (the first row), En null mutant (the second row), Hh and Ci null mutants (the third row), and the TGF-β and Wg null mutants (the fourth row). The pink lines atradii0.3 and 0.6 are used to separate the radius into inner ring between 0 and 0.3, middle ring between 0.3 and 0.6, and outer ring between 0.6 and 1. The expected pigment in the corresponding region is shown at the top of each region. (C) displays the species with the corresponding eyespot pattern marked in the yellow boxes. The first row shows the eyespot with white inner ring, black middle ring, and yellow outer ring in Bicyclus anynana. The second row provides the eyespot with black inner ring in Vanessa ataianta.

The third row includes two cases for the white inner ring: one in Vanessa atalanta and the other in Vanessa altissima. The fourth row shows the white inner ring and black middle ring in Chlosyne nycteis. A cartoon eyespot pattern representing each phenotypic case is shown above each arrow between (B) and (C). The experiment and prediction listed below each arrow between (B) and (C) represent that the expected eyespot pattern is with and without experimental evidence, respectively.

**,QQHUULQJ** **0LGGOHULQJ** **2XWHUULQJ**

'OO (Q

+K &L 7*)ߚ (Q

&L :J

**5DGLXV**_{}

Figure 2. System network of the model. Initially, for the cells in the inner ring, the protein Dll triggers the gene expression of En to generate the mRNA and protein of En [40]. The generated En protein then triggers the gene expression of Hh [23], and then the produced Hh protein inhibits the gene expression of Ci in the inner ring [20, 67]. Next, the Hh protein diffuses to cell member to bind with the receptor patched on the cells in the middle ring to activate the phosphorylation of Ci protein in the middle ring [23]. In the middle ring, the phosphated Ci protein triggers the gene expression of TGF-β [23], and then the generated TGF-β protein diffuses to the outer ring to initiate the

autoregulation of En protein [23]. Meanwhile, the TGF-β protein diffuses to the inner ring to inhibit the production of Hh protein [10]resulting in promoting Ci protein in the inner ring. Next, in the inner ring, the generated Ci protein binds with Wg protein to trigger the Wg signaling pathway [11].

The triggered Wg signaling pathway diffuses to the outer ring [55]to maintain the expression of En activated by the TGF-β [29, 51].

The equation for En protein (E)

188

The equation for the concentration of En protein is described by

189

∂E

∂t = αE· χin(r)

| {z }

source from Dll

+ (λET · χout(r)) · (W · χout(r))

| {z }

interaction between Wg and TGF-β

− µEE

| {z }

degradation

. (5)

The first term in the right-hand side of Eq. (5) represents a constant source of En protein from Dll in

190

the inner ring [40]. The second term accounts for the autoregulation of En protein [23]maintained by

191

the interaction between TGF-β protein and Wg protein in the outer ring [29, 51]. The last term is the

192

degradation of En protein.

193

The equation for Hh protein (H)

194

The concentration of Hh protein satisfies the following equation

195

∂H

∂t = dH

1
r^{2}

∂

∂r(r^{2}∂H

∂r)

| {z }

diffusion

+ (αHE · χin(r))

| {z }

promotion triggered by En

/(1 + kHT · χin(r))

| {z }

inhibition by TGF-β

− µHH

| {z }

degradation

. (6)

The first and last terms account for the diffusion and degradation of Hh protein. The second term shows

196

that the Hh protein is triggered by En protein [23]and is inhibited by TGF-β protein [10]in the inner

197

ring.

198

Table 2. Parameters of the whole model.

Parameter Description Value Unit Reference

dH diffusion rate of Hh 2.97017 × 10^{−7} cm^{2}/min [12] & Table 6

dT diffusion rate of TGF-β 3 × 10^{−7} cm^{2}/min [12] & Table 6

dW diffusion rate of Wg 2.91519 × 10^{−7} cm^{2}/min [12] & Table 6

µ_{E} degradation rate of En 3.85082 × 10^{−4} /min [35] & Table 7

µ_{H} degradation rate of Hh 1.38629 × 10^{−2} /min [18, 31, 56] & Table 6

µ_{C} degradation rate of Ci 9.24196 × 10^{−3} /min [4] & Table 7

µ_{T} degradation rate of TGF-β 5.77623 × 10^{−3} /min [60] & Table 6

µW degradation rate of Wg 1.15525 × 10^{−3} /min [73] & Table 7

N0 amount of Wg protein 1 kD/cm estimated

R0 radius of eyespot 0.094 cm [55]

α_{E} production rate of En in inner ring (3.56201 × 10^{−2})N_{0} kD/cm [55]& estimated
α_{H} production rate of Hh in inner ring 2.18247 × 10^{−2} /min [55]& estimated
α_{C} production rate of Ci in middle ring 1.20978 × 10^{−3} /min [55]& estimated
α_{T} production rate of TGF-β in middle ring 1.89385 × 10^{−5} /min [55]& estimated
α_{W} production rate of Wg in inner ring 9.46923 × 10^{−5} /min [55]& estimated

λE production rate of En in outer ring ^{3.56201×10}_{N} ^{−2}

0 cm/kD/min [55]& estimated λC production rate of Ci in inner ring 352.35N0 kD/cm [55]& estimated

k_{H} half-saturation of Hh 1/(4N_{0}) cm/kD estimated

k_{C} half-saturation of Ci 1/(932N_{0}) cm/kD estimated

The equation for Ci protein (C)

199

The equation for the concentration of Ci protein is given by

200

∂C

∂t = λC

|{z}

production of Ci

· χin(r)
1 + k_{C}H

| {z }

inhibition by Hh in inner ring

+ αC

|{z}

production of Ci

· H · χmid(r)

| {z }

promotion by Hh in middle ring

− µCC

| {z }

degradation

. (7)

The first term represents the production of Ci protein inhibited by Hh protein in the inner ring [20, 67].

201

The second term shows the production of Ci protein trigged by Hh protein in the middle ring [23].

202

The equation for TGF-β protein (T )

203

We model the dynamics of TGF-β protein by the equation

204

∂T

∂t = d_{T} 1
r^{2}

∂

∂r(r^{2}∂T

∂r)

| {z }

diffusion

+ α_{T}C · χ_{mid}(r)

| {z }

activation of TGF-β in middle ring

− µ_{T}T

| {z }

degradation

. (8)

The first and the last terms account for the diffusion and degradation of TGF-β protein. The second

205

term represents the activation of TGF-β protein by Ci protein in the middle ring [23].

206

The equation for Wg protein (W )

207

The concentration of Wg protein satisfies the following equation

208

∂W

∂t = dW

1
r^{2}

∂

∂r(r^{2}∂W

∂r)

| {z }

diffusion

+ αWC · χin(r)

| {z }

activation by Ci

− µWW

| {z }

degradation

. (9)

The first and last terms account for the diffusion and degradation of Wg protein. The second term shows

209

the activation of Wg protein by Ci protein in the inner ring [11].

210

The detailed molecular information to support Eqs. (5)-(9) comes from experimental data of insects

211

with wings (mostly from Drosophila) [10, 11, 20, 23, 29, 40, 51, 67]. Since insects with wings, Drosophila,

212

butterflies, and moths havea shared ancestor that also had wings [2], they share numerous similarities of

213

wing structure and morphology, wing development, and molecular processes, including genetic architec-

214

ture and mechanisms of gene regulation [2]. Thus, we transfer the molecular knowledge fromwell-studied

215

insects with wings to butterflies to estimate the parameter values for butterflies and then use sensitivity

216

analysis to study how these parameter values correlate to the eyespot pattern.

217

Initial condition (IC).

218

Initially, only En is presented in the inner ring due to the Dll gene expression. Therefore, we assume

219

that En in the inner ring is a decreasing function with respect to the radius r with the maximum

220

occurring at the center, i.e., r = 0, and En maintains the minimum in the middle and outer rings. The

221

initial conditions for other variables are given based on their steady states shown in the Appendix. So

222

we have

223

E(r, 0) = 185 × 10^{−5}+ 10^{−2} 1 −10
3r

!

χ_{[0,0.3]}(r)
H(r, 0) = 233 × 10^{−5}

C(r, 0) = 61 × 10^{−5}
T (r, 0) = 10^{−5}
W (r, 0) = 10^{−5},

(10)

for 0 ≤ r ≤ 1.

224

Boundary condition (BC) - Neumann BC.

225

All dependent variables are radially symmetric and have no flux at the boundary of the eyespot. In

226

other words, we have

227

∂η

∂r(1, t) = 0 (11)

for η = {E, H, C, T, W } and t ≥ 0. On the other hand, to guarantee that the solution is regular, i.e., lim sup

r→0^{+}

1 r

∂u

∂r(r, t) < ∞, it is necessary that

228

∂η

∂r(0, t) = 0 (12)

for t ≥ 0.

229

Combining above equations, we obtain the following system of PDEs

230

∂E

∂t = αE· χin(r) + (λET · χout(r)) · (W · χout(r)) − µEE

∂H

∂t = dH

1
r^{2}

∂

∂r(r^{2}∂H

∂r) + (αHE · χin(r))/(1 + kHT · χin(r)) − µHH

∂C

∂t = λCχin(r)

1 + kCH + αCH · χmid(r) − µCC

∂T

∂t = dT

1
r^{2}

∂

∂r(r^{2}∂T

∂r) + αTC · χmid(r) − µTT

∂W

∂t = dW

1
r^{2}

∂

∂r(r^{2}∂W

∂r) + αWC · χin(r) − µWW

(13)

with initial condition (10) and homogeneous Neumann boundary condition (11) and (12). For the model

231

(13), as mentioned in Fig. 2, there are experimental data from insects with wings to support the qual-

232

itative interactions among the components [10, 11, 20, 23, 29, 40, 51, 67]. Thus, the assumptions of the

233

model (13) are: (i) the symmetric shape of the eyespot, initial and boundary conditions, while all of these

234

assumptions are close to the real situation; (ii) the similarity of the structure development and molec-

235

ular processes between insects with wings and butterflies; and (iii) the estimated parameter values of

236

k_{H}, k_{C}, µ_{E}, µ_{H}, and N_{0}. Notice that different value of N_{0}affects the absolute values of α_{E}, λ_{E}, λ_{C}, k_{H},

237

and k_{C}, but the relative values amount these five values is fixed. Since the assumptions (ii) and (iii) affect

238

the values of k_{H}, k_{C}, µ_{E}, µ_{H}, and N_{0}, we will perform sensitivity analysis to demonstrate that these

239

five parameter values do not have significant effect on the model outcome. Thus, when the values of these

240

five parameters in butterflies become available, the simulation outcome will not be changed or we will be

241

able to predict the simulation outcome easily.

242

### Numerical simulation

243

In this section, we use the model (13) with conditions (10), (11), (12) to simulate the patterns of eyespots

244

in wild type(Bicyclus anynana) and null mutants of butterflies. We use the forward Euler method with

245

time step dt = 0.01 minute and dr = 0.001/R_{0} and run the program in MATLAB. We want to remark

246

that due to (12) we drop the term 1 r

∂η

∂rand take

247

1
r^{2}

∂

∂r(r^{2}∂η

∂r) ≈∂^{2}η

∂r^{2}

in our simulation. Since the dropped term is of lower order, the qualitative behaviors of the solutions

248

will not be affected. In the following, the notations En, Hh, Ci, TGF-β, and Wg account for En, Hh, Ci,

249

TGF-β, and Wg proteins, respectively.

250

Model validation by using the wild type

251

In this subsection, we first validate the mathematical model (13) by comparing the numerical simulation

252

results with the patterns of eyespots in wild type species. Notice that the high concentration of En

253

triggers the generation of white and yellow pigments in the inner and outer rings, respectively. A low

254

concentration of En, the high concentration of Ci activates the production of black pigment. However,

255

when both of En and Ci have high concentrations at the same location, the cells only generate white

256

pigments. Additionally, since the Dll protein triggers the En expression in the inner ring, it takes around

257

16 hours to generate all pigments in the corresponding rings. To mimic the wild type eyespot pattern

258

with white, black, and yellow pigments in the inner, middle, and outer rings, we expect that, at 16 hours,

259

(i) both of En and Ci have high concentrations in the inner ring for white pigment; (ii) only Ci has a

260

high concentration in the middle ring for black pigment; and (iii) only En has a high concentration in

261

the outer ring for yellow pigment.

262

Fig. 3 shows the time series of the simulation results of model (13) for wild type at time {6 min, 30 min, 1 hr, 2 hr, 4 hr, 6 hr, 8 hr, 16 hr}.

First, the initial peak of En in the inner ring induces the peak of Hh in the inner ring during the early

263

period (around 1 hour). The induced Hh immediately causes an increasein the concentration of Ci in

264

the middle ring and inhibits the concentration of Ci in the inner ring, resulting ina singlepeak of Ci in

265

the middle ring (around 1 hour). The Cigeneratedin the middle ringthenproduces the peak of TGF-β

266

in the middle ring (at about 4 hours). The diffusion of TGF-β immediately increases the production of

267

Ci in the inner ring to generate the peaks of Ci and Wg in the inner ring (around 4 hours). Finally,

268

both TGF-β and Wg diffuse to the outer wing to generate the peak of En there (around 4 hours). All

269

theseinteraction are maintainedsuch that the amplitudes of these peaks keepincreasing. Eventually, En

270

develops two distinctpeaks in the inner and outer rings, Hh has a peak in the inner ring, Ci has peaks

271

in the inner and middle rings, TGF-β has a peak in the middle ring, and Wg has a peak in the inner

272

ring. The high concentrations of En and Ci in the inner ring trigger cells to produce white pigment in the

273

inner ring. The high concentrations of Ci in the middle ring and En in the outer ring activate precursor

274

cells to generate black pigment in the middle ring and yellow pigment in the outer ring, respectively.

275

Additionally, experimental observations show that (i) the concentration of En is high in the inner and

276

outer rings [9], (ii) the concentration of Hh is high in the inner ring [23], (iii) the concentration of Ci is

277

high in the inner and middle ring [23], (iv) the concentration of TGF-β is high in the middle ring [41], and

278

(v) the concentration of Wg is high in the inner ring [41], at the end of eyespot formation process, i.e., 16

279

hours. Therefore, our mathematical model (13) generates the concentration profiles of these five genes in

280

wild type, which is in accordance with the experimental observation at 16 hours. However, experimental

281

data of the temporal dynamics of these five gene expression concentrations over time are unavailable so

282

far. Thus, we are unable to compare our temporal dynamics simulation (Fig. 3A-G) with experimental

283

observation.

284

Numerical predictions for null mutants

285

In this subsection, we use the model (13) to numericallyinvestigate the eyespot pattern when the com-

286

ponents En, Hh, Ci, TGF-β, and Wg are knockout separately, i.e., the null mutants, to motivate future

287

experiments for validation.

288

First, we consider the En null mutants (namely, knockout the En in the model (13)) that we set

289

E(r, t) ≡ 0, for 0 ≤ r ≤ 1 and for all t ≥ 0. The simulation result of the En null mutants is shown in Fig.

290

4. In Fig. 4, when the En is knockout, the concentration of Hh is null everywhere, leading to the result

291

that there is no Ci in the middle ring to generate the black pigment and only a low concentration of

292

TGF-β appears in the middle ring (Notice that the maximal value of TGF-β is around 2 × 10^{−9}kD/cm

293

which is neglectable, comparing to the wild type in Fig. 3H. The lack of TGF-β eliminates the peak of

294

En in the outer ring, even though the profile of Wg is similar to the wild type. The peaks of En in the

295

inner and outer rings disappear, but Ci still has a peak in the inner ring. Thus, the En null mutants

296

only generate a single inner ring with black pigment, since Ci generates black pigment under the absence

297

of En. This kind of degenerated eyespot is observed from the butterfly, Vanessa atalanta. The Vanessa

298

atalanta has hindwing eyespots where only the inner ring with black pigment is present as shown in Fig.

299

9A. Hence, we hypothesize that En deficiency shifts the black pigment to the inner ring and then causes

300

the degenerate eyespot pattern with a single black spot.

301

Fig. 5 shows the simulation result of Hh null mutants by setting H(r, t) ≡ 0, for 0 ≤ r ≤ 1 and t ≥ 0.

302

The En is normally expressed by Dll in the inner ring, but lack of Hh blocks the interaction between En

303

and Ci. Thus, the profiles of Ci, TGF-β, and Wg are the same as the ones shown in En null mutants,

304

i.e., Fig. 4, that the peak of Ci in the middle ring and the peak of En in the outer ring vanish. Hence,

305

the Hh null mutants only have one peak of En and one peak of Ci in the inner ring, resulting in a single

306

white spot.

307

For the Ci null mutants, we set C(r, t) ≡ 0, for 0 ≤ r ≤ 1 and t ≥ 0. The simulation result is displayed

308

in Fig. 6. Since the En works normally, the peaks of En and Hh in the inner ring exist. However, the Ci

309

is vanished everywhere, so there is no peaks of Ci and none of TGF-β and Wg are generated resulting

310

in the loss of En peak in the outer ring. Hence, the Ci null mutants only generate one peak of En in the

311

inner ring resulting in a single white spot.

312

Based on the results in Figs. 5 and 6, lack of Hh and lack of Ci generate the same degenerated eyespot

313

pattern: a single white spot, which is observed from the species: Vanessa atalanta (see Fig. 9A) and

314

Vanessa altissima (see Fig. 9B). Thus, we hypothesize that the degenerated single white spot pattern in

315

these two species is caused by the deficiency of Hh or Ci.

316

Figure 3. Time series of simulation results for wild type eyespot pattern. (A)-(H) show the simulation results of the model (13) at {6 min, 30 min, 1 hr, 2 hr, 4 hr, 6 hr, 8 hr, 16 hr},

respectively. In each figure, the first, second, third, fourth, and fifth rows display the concentrations of
En, Hh, Ci, TGF-β, and Wg, respectively. The horizontal and vertical axes represent the radius with
unit R0(= 0.094 cm) and the concentration of protein with unit kD/cm. For the radius, regions
[0, 0.3], [0.3, 0.6], and [0.6, 1] represent the inner ring Ωin, middle ring Ωmid, and outer ring Ωout,
respectively. (H) shows the final stage of the eyespot formation, i.e., at 16 hours, that (i) the maximal
values of En in Ω_{in} and Ω_{out} are at 3.64603 × 10^{−2} kD/cm and 2.96234 × 10^{−2} kD/cm; (ii) the
maximal values of Ci in Ω_{in}and Ω_{mid} are at 3.33815 × 10^{−1} kD/cm and 2.69242 × 10^{−1} kD/cm; and
(iii) the maximal values of Hh in Ω_{in}, TGF-β in Ω_{mid}, and Wg in Ω_{in}are at 5.04047 × 10^{−1} kD/cm,
1.23412 × 10^{−5} kD/cm, and 1.49298 × 10^{−3} kD/cm, respectively.

Next, in Fig. 7, we consider the TGF-β null mutants with T (r, t) ≡ 0, for 0 ≤ r ≤ 1 and t ≥ 0. The

317

profiles of En, Hh, and Ci in the inner and middle rings are similar to the wild type. However, knockout

318

of TGF-β blocks the production of Wg such that no interaction between TGF-β and Wg in the outer

319

ring resulting in the loss of En peak in the outer ring. Therefore, the TGF-β null mutants have peaks of

320

En and Ci in the inner ring to produce white pigments in the inner ring and one peak of Ci in the middle

321

ring to generate black pigments in the middle ring. However, cells lose the ability to generate the yellow

322

Figure 4. Simulation of En null mutants. (A)-(E) show the concentrations of En, Hh, Ci, TGF-β,
and Wg, respectively, at 16 hours,inEn null mutants. The horizontal and vertical axes represent the
radius with unitR0(= 0.094 cm)and concentration of protein with unit kD/cm. The maximal value of
TGF-β is at 2.41413 × 10^{−9} kD/cm.

pigments in the outer ring due to the loss of En peak in the outer ring.

323

Finally, we consider the Wg null mutants with W (r, t) ≡ 0, for 0 ≤ r ≤ 1 and t ≥ 0 in Fig. 8. The

324

profiles of En, Hh, Ci, and TGF-β in the inner and middle rings are similar to the wild type. However,

325

knockout of Wg loses the interaction between TGF-β and Wg in the outer ring, such that the En peak

326

in the outer ring disappears resulting in loss of the yellow pigments. Hence, the Wg null mutants have

327

the same eyespot pattern as the TGF-β null mutants that the outer yellow ring disappears. This type of

328

degenerated eyespot pattern, losing the outer yellow ring, can be observed from the butterfly, Chlosyne

329

nycteis (see Fig. 9C). Combining the results from Figs. 7 and 8, we conjecture that the degenerated

330

pattern with losing the outer yellow ring is caused by loss of TGF-β or Wg signaling.

331

From the above simulations,we predict the following three types of degenerated patterns in knockout

332

mutantsand the results are summarizedin Table 3:

333

(i) deficiency of En causes a single black spot, which can be observed from Vanessa atalanta;

334

(ii) deficiency of Hh or Ci generates a single white spot, which can be observed from Vanessa atalanta

335

and Vanessa altissima;

336

(iii) deficiency of TGF-β or Wg loses the outer yellow ring, which can be observed from Chlosyne nycteis.

337

Additionally, the temporal dynamics of these five knockout null mutations are similar to the wild type

338

case shown in Fig. 3 that the time series of Figs. 4-8 show similar profiles of these five components during

339

the whole process. This means that the stable patterns appear at the beginning and maintain during the

340

whole process, so there is no bifurcation or different profiles appear in these five mutations.

341

Currently, there are a lack of knockout experiments to validate our numerical knockout predictions, so

342

future experiments are required for validation. There are two possible approaches for future experiments.

343

The first approach is to study the phenotype of these null mutants by using the genotype, which includes

344

two types of knockout experiments. One is to knockout or severely knock down one of the components

345

Figure 5. Simulation of Hh null mutants. (A)-(E) show the concentrations of En, Hh, Ci, TGF-β,
and Wg, respectively, at 16 hours,inHh null mutants. The horizontal and vertical axes represent the
radius with unitR0(= 0.094 cm)and concentration of protein with unit kD/cm. The maximal value of
TGF-β is at 2.41413 × 10^{−9} kD/cm.

{En, Hh, Ci, T GF − β, W g} in the embryo [30], through the use of CRISPR, morpholinos, RNAi, or

346

dominant negative viral vector constructs. However, these genetic components play essential roles in the

347

early stages of butterfly development [30], so removing any of these components will lead to embryonic or

348

larval lethality. Thus, it is difficulty to collect gene expression data in late larval and pupal butterfly wings

349

discs form this type of knockout experiment. However, this problem could be solved if the gene knockout

350

process can be performed later in development, perhaps in the late fourth instar or early fifth instar larva

351

stage immediately prior to eyespot determination. The other type of knockout experiment is to knockout

352

essential enzymes for pigment synthesis [30, 71] such that it will not affect essential tissue and organ

353

formation but cells lose the ability to generate the pigments in the corresponding rings. However, this

354

kind of mutant is different from our simulation setting because none of the components of eyespot ring

355

specification are eliminated, and hence it cannot be used to validate our prediction results. An alternative

356

to study of knockout or knock down experiments is to study the genotype by using the phenotype from

357

the selected lines. For the existing experiments in Bicyclus anynana, there are two selected lines: one

358

is no black ring (c.f. [9]) and the other is no yellow ring (c.f. [3] and Chlosyne nycteis in Fig. 9D).

359

The difference of gene expression between the selected lines and wild type can be used to study the key

360

factors for generating the deficient eyespot pattern.

361

### Sensitivity analysis

362

(We change everything in this Section.)

363

In this section, we perform the sensitivity analysis created by S. Marino et al. in [36] to investigate

364

the robustness of the model outcomes and prevent the over fitting issue, by analyzing how the parameter

365

values affect the wild type (Bicyclus anynana) eyespot pattern (namely, the peaks of En and Ci in

366

different rings) and the variations of all components. We will first apply the sensitivity analysis on three

367

Figure 6. Simulation of Ci null mutants. (A)-(E) show the concentrations of En, Hh, Ci, TGF-β,
and Wg, respectively, at 16 hours,inCi null mutants. The horizontal and vertical axes represent the
radius with unitR_{0}(= 0.094 cm)and concentration of protein with unit kD/cm.

Figure 7. Simulation of TGF-β null mutants. (A)-(E) show the concentrations of En, Hh, Ci,
TGF-β, and Wg, respectively, at 16 hours,in TGF-β null mutants. The horizontal and vertical axes
represent the radius with unitR_{0}(= 0.094 cm)and concentration of protein with unit kD/cm.

Figure 8. Simulation of Wg null mutants. (A)-(E) show the concentrations of En, Hh, Ci, TGF-β, and Wg, respectively, at 16 hours,inWg null mutants. The horizontal and vertical axes represent the radius with unitR0(= 0.094 cm)and concentration of protein with unit kD/cm.

Table 3. Eyespot patterns for different null mutants. In the species row, A-D represent the wild type of Bicyclus anynana, Vanessa atalanta, Vanessa altissima, and Chlosyne nycteis, respectively. The named pattern row shows the cartoon of the expected eyespot pattern. The actual eyespot patterns are shown in Fig. 9.

Mutant type wild type null En null Hh null Ci null TGF-β null Wg
En peak in Ω_{in} presence absence presence presence presence presence

pigment in Ωin white black white white white white

Ci peak in Ω_{mid} presence absence absence absence presence presence

pigment in Ωmid black null null null black black

En peak in Ωout presence absence absence absence absence absence

pigment in Ωout yellow null null null null null

pattern

species A B B, C B, C D D

cases to study the wild type pattern: (i) the peak of En in the outer ring, (ii) the peak of Ci in the

368

inner ring, and (iii) the peak of Ci in the middle ring. Next, for the model robustness, we perform the

369

sensitivity analysis on all components in all rings to broadly investigate how the parameter values affect

370

the gene expression pattern.

371

The concept of sensitivity analysis mentioned in [36] is to evaluate how the uncertainty and variations

372

in model outputs are correlated to parameter values, by using the Latin hypercube sampling (LHS)

373

and partial rank correlation coefficient (PRCC). For each parameter, the LHS is a sampling method

374

Figure 9. Images of butterfly. (A)-(D) show the images of butterflies: wild type of Bicyclus anynana, Vanessa atalanta, Vanessa altissima, and Chlosyne nycteis, respectively. In (A), the yellow box indicates the eyespot of wild type. In (B), the yellow and pink boxes show the black inner ring and white inner rings, respectively, in Vanessa atalanta. The yellow box in (C) display the white inner ring in Vanessa altissima, and the pink box in (D) indicates the white inner ring and black middle ring in Chlosyne nycteis.

that generates uniform parameter value distributions divided into N equal probability intervals, where

375

N is the sample size. Each interested parameter will be sampled independently by using LHS. All the

376

samples are collected to generate a set {P1, P2, · · · , PN} and each Piincludes the values for all interested

377

parameters. Next, substitute each set Pi into the parameter values of the model to generate the model

378

outcomes {y1, y2, · · · , yN}. The value of PRCC between the parameter values {P1, P2, · · · , PN} and model

379

outcomes {y1, y2, · · · , yN} shows the robust sensitivity for their nonlinear and monotonic relationships.

380

Thus, a parameter with positive PRCC to the model outcome and p-value smaller than 0.05 represents

381

that the model outcome increases as the value of the parameter increases, whereas a parameter with

382

negative PRCC to the model outcomes and p-value smaller than 0.05 accounts for an opposite result that

383

the model outcome decreases as the value of the parameter increases. However, for a parameter with

384

small |P RCC| and p-value larger than 0.05, then the parameter value does not have significant effect on

385

the model outcome. Thus, the sensitivity analysis can also be used to study the over fitting issue. If

386

most parameters are with small |P RCC| and p-value larger than 0.05, then the model outcome is not

387

sensitive to most parameter values which means that the dynamics of the model is robustness.

388

To analyze how the parameter values affect the peak appearance in each ring, we define the following

389

functions

390

Xˆin(t) := max{X(r, t) : 0 ≤ r ≤ 0.3},
Xˆ_{mid}(t) := max{X(r, t) : 0.3 ≤ r ≤ 0.6},

Xˆout(t) := max{X(r, t) : 0.6 ≤ r ≤ 1}, and

391

X_{in}(t) =

Xˆin(t)

X(0.3, t), X_{mid}(t) =

Xˆmid(t)

X(0.3, t) + X(0.6, t), X_{out}(t) =

Xˆout(t)

X(0.6, t), (14) with X ∈ {E, H, C, T, W }. Since we only focus on the end of the development, we take t = 16 hours

392

in this section. By using the functions in Eq. (14), if Xj(t) > 1 with X ∈ {E, H, C, T, W } and

393

j ∈ {in, mid, out}, then there exists at least one peak of X in the ring j, at time t. For each parameter,

394

we generated 10000 samples individually, via the Latin hypercube sampling, with dt = 0.01 minutes and

395

dr = 0.001 × R0. We choose these parameters in the range from 0.5 to 2 fold of their baseline values.

396

The Table 4 shows the baselines, ranges, and units of the parameters for each parameter.

397

Our simulations show that the following functions are always smaller than one under all parameter samples,

E_{mid}(t), H_{mid}(t), H_{out}(t), C_{out}(t), T_{in}(t), T_{mid}(t), T_{out}(t), W_{mid}(t), W_{out}(t).

This result indices that there is no peak of Hh, Ci, TGF-β, and Wg in the outer ring and TGF-β in

398

the inner ring, and we cannot make any conclusion of the peak in the middle ring for En, Hh, TGF-β,

399

and Wg. Therefore, in the following, we only focus on the En peaks in the inner and outer rings, Hh

400

peak in the inner ring, Ci peaks in the inner and middle rings, and Wg peak in the inner ring. The

401

PRCCs and the p-values of parameters for these cases are shown in Table 5. For the considered Xj(t), a

402

parameter with a negative (resp. positive) PRCC and p-value smaller than 0.05 indicates that increasing

403

this parameter will decrease (resp. increase) the ratio Xj(t) and hence reduces (resp. increases) the

404

chance to generate the peak of X in the ring j.

405

In the following, we will use the PRCCs in Table 5 to investigate i) how the parameter values affect

406

the peaks in wild type, ii) how the diffusion and degradation rates of each gene affect the peaks, and iii)

407

how the parameter values affect the remainder peaks.

408

Wild type pattern.

409

The eyespot pattern in the wild type requires a peak of En and a peak of Ci in the inner ring for white

410

pigment, a peak of Ci in the middle ring for black pigment, and a peak of En in the outer ring for yellow

411

pigment. According to the initial condition, the peak of En in the inner ring always exists. Thus, we only

412

need to consider three cases: (i) the peak of En in the outer ring, (ii) the peak of Ci in the inner ring,

413

and (iii) the peak of Ci in the middle ring. We choose the parameters λE and αE for case (i), λC, αH,

414

and αC for case (ii), and αC, µE, αE, αH, λC, and dH for case (iii). We then use the PRCCs with

415

p-value corresponding to the ratios Eout(t) in case (i), Cin(t) in case (ii), and Cmid(t) in case (iii), at

416

t = 16 hours.

417

In case (i), the PRCC of αE and λE are positively correlated to Eout(t)that increasing αE or λE

418

promotes the generation of En peak in the outer ring. Increasing αE generates more En as a source

419

of other components and hence increases the amount of Wg to generate more En in the outer ring.

420

Increasing λE enhances the effect from the interaction between TGF-β and Wg in the outer ring, such

421

that En is more sensitive to TGF-β and Wg for generating the peak in the outer ring.

422

In case (ii), the PRCCs of αHand αCare negatively correlated and λCis positive correlated to Cin(t).

423

Increasing α_{H} promotes the production of Hh in the inner ring resulting in more inhibition on Ci in the

424

inner ring. On the other hand, increasing α_{C}enhances the production of Ci in the middle ring. However,

425

the negative feedback of Ci then inhibits the production of Hh in the inner ring and hence increases the

426