Mathematical modeling of the eyespots in butterfly wings
1
Kang-Ling Liao1,2,∗ Wei-Chen Chang3, Jeffrey M. Marcus2, Jenn-Nan Wang4
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 mm2 [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 + e30(r−0.15), (2)
χmid(r) = 1.5 × 1
1 + e−30(r−0.4)× 1
1 + e30(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 r2
∂
∂r(r2∂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 cm2/min [12] & Table 6
dT diffusion rate of TGF-β 3 × 10−7 cm2/min [12] & Table 6
dW diffusion rate of Wg 2.91519 × 10−7 cm2/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)N0 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×10N −2
0 cm/kD/min [55]& estimated λC production rate of Ci in inner ring 352.35N0 kD/cm [55]& estimated
kH half-saturation of Hh 1/(4N0) cm/kD estimated
kC half-saturation of Ci 1/(932N0) 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 + kCH
| {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 = dT 1 r2
∂
∂r(r2∂T
∂r)
| {z }
diffusion
+ αTC · χmid(r)
| {z }
activation of TGF-β in middle ring
− µTT
| {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 r2
∂
∂r(r2∂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 r2
∂
∂r(r2∂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 r2
∂
∂r(r2∂T
∂r) + αTC · χmid(r) − µTT
∂W
∂t = dW
1 r2
∂
∂r(r2∂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
kH, kC, µE, µH, and N0. Notice that different value of N0affects the absolute values of αE, λE, λC, kH,
237
and kC, but the relative values amount these five values is fixed. Since the assumptions (ii) and (iii) affect
238
the values of kH, kC, µE, µH, and N0, 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/R0 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 r2
∂
∂r(r2∂η
∂r) ≈∂2η
∂r2
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−9kD/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 Ωinand Ω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 Ωinare 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 unitR0(= 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 unitR0(= 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
Xin(t) =
Xˆin(t)
X(0.3, t), Xmid(t) =
Xˆmid(t)
X(0.3, t) + X(0.6, t), Xout(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,
Emid(t), Hmid(t), Hout(t), Cout(t), Tin(t), Tmid(t), Tout(t), Wmid(t), Wout(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 αCenhances 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
