Mathematical models for the signaling pathway of the G-protein coupled receptor EP2

Ti le: Ma hema ical models o he signaling pa hway o he G-p o ein
coupled ecep o EP2
Au ho : Paula Gómez López
Ad iso : Gemma Hugue Casades
Co-ad iso : S e anie Sonne
Depa men : Depa men o Ma hema ics
Academic yea : 2024/2025
Mas e o Science in
Ad anced Ma hema ics and
Ma hema ical Enginee ing
Uni e si a Poli `ecnica de Ca alunya
Facul a de Ma em`a iques i Es ad´ıs ica
Mas e in Ad anced Ma hema ics and Ma hema ical Enginee ing
Mas e ’s hesis
Ma hema ical models o
he signaling pa hway o he
G-p o ein coupled ecep o EP2
Paula G´omez L´opez
Supe ised by D . Gemma Hugue Casades and D . S e anie Sonne
Janua y, 2025
This mas e ’s hesis was ca ied ou wi h he suppo o he E asmus+ In e nship p og am om Sep em-
be o Decembe 2024 a Radboud Uni e si y, unde he supe ision o D . S e anie Sonne .

Fi s o all, I would like o hank S e anie Sonne o guiding me h ough his hesis and o answe ing
all my ques ions. I deeply app ecia e Ma iya P ashnyk o all he online mee ings and he help wi h he
code. Thanks o my pa en s o always being he e, o Daphne o all he suppo , and o my iends—bo h
old and new— o making his expe ience so much be e .
Abs ac
This hesis de elops ma hema ical models o s udy he EP2 signaling pa hway, a ecep o belonging o he
G-p o ein-coupled ecep o (GPCR) amily in ol ed in cellula communica ion. By combining ma hema ical
echniques and biological insigh s, we aim o be e unde s and he dynamics o EP2 signaling and explo e
whe he spa ial e ec s play an impo an ole.
Key con ibu ions include he de elopmen o an ODE model o desc ibe he dynamics o cAMP p o-
duc ion in EP2 signaling and i s ex ension o inco po a e ligand- ecep o in e ac ions. Spa ial e ec s a e
subsequen ly in es iga ed h ough eac ion-di usion models. Finally, a no el class o ligand- ecep o -based
Tu ing models is in oduced o explo e ecep o clus e ing ia spa ial pa e n o ma ion, p o iding a heo-
e ical amewo k applicable o a ious ecep o sys ems.
This wo k highligh s he powe o ma hema ical modeling in biological signaling and o e s ools ha
ex end beyond he EP2 ecep o , in e connec ing ma hema ics and biology.
Keywo ds
Ma hema ical modeling, G-p o ein coupled ecep o (GPCR), EP2 signaling, ligand- ecep o dynamics,
Tu ing pa e ns, Schnakenbe g ype kine ics.
1
Con en s
1 In oduc ion 3
1.1 Biological backg ound ..................................... 4
1.2 Expe imen al da a ....................................... 6
2 Modeling he signaling pa hway o EP2 8
2.1 G-p o ein cycle ......................................... 8
2.2 cAMP p oduc ion ....................................... 12
2.3 ODE model o cAMP p oduc ion in EP2 ecep o signaling ................. 14
2.3.1 Simula ions ...................................... 14
2.4 Modeling ligand- ecep o dynamics .............................. 16
2.4.1 Simula ions ...................................... 18
3 Spa ially he e ogeneous model o EP2 signaling 20
3.1 1D eac ion-di usion model in a c oss-sec ion o he cell .................. 21
3.1.1 Simula ions ...................................... 22
3.2 1D eac ion-di usion model on he cell memb ane ...................... 27
3.2.1 Simula ions ...................................... 29
4 Ligand- ecep o based Tu ing models 32
4.1 Gene al condi ions o di usion-d i en ins abili y ....................... 33
4.1.1 Linea s abili y in he absence o di usion ...................... 33
4.1.2 Di usion d i en ins abili y .............................. 34
4.2 Ligand ecep o model ..................................... 38
4.3 Ligand- ecep o based Tu ing models ............................. 40
4.3.1 Gene alized Schnakenbe g model ........................... 42
4.3.2 Model wi h no eedback ecep o p oduc ion ..................... 51
5 Conclusion and u u e wo k 54
6 Re e ences 55
A Codes 57
2
Mo eo e , ollowing [12], we include he e e sible Reac ion 1in ou model, whe eas [10] ea ed his
eac ion as one-way. The o he eac ions could also be modeled as e e sible, depending on he molecula
concen a ions and eac ion a es.
In [10], he Reac ions 2,3a e ea ed as enzyma ic, ha is, he ac i a ed ecep o ca alyzes he dis-
socia ion o he G-p o ein, while GAP ca alyzes GTP hyd olysis, modeled by Michaelis-Men en kine ics.
In con as , we assume ha Reac ions 1and 3a e non-enzyma ic and only Reac ion 2 ollows Michaelis-
Men en kine ics. Reac ions 1and 3a e modeled using he Law o Mass Ac ion, which s a es ha he a e
o a eac ion is p opo ional o he p oduc o he concen a ions o he eac an s [18].
Le [·] deno e he concen a ion o a chemical subs ance. Mo eo e , GαGDP
sβγ is deno ed by αsβγ,
GαGDP
sby Gαsand GαGTP
sby α∗
s. Then, he eac ion a es a e de ined as ollows:
1. Associa ion/dissocia ion o Gαsβγ: The e e sible eac ion, ep esen ing he binding and unbinding
o GαGDP
sand Gβγ, is no enzyme-ca alyzed and is gi en by:
V1=k3[αs][βγ]−k4[αsβγ],
whe e k3and k4a e he Gαsβγ associa ion and dissocia ion a e espec i ely.
2. Ac i a ion o Gαs: The enzyma ic eac ion, whe e he ac i a ed ecep o , [EP2∗], ca alyzes he
dissocia ion o Gαsβγ in o GαGTP
sand Gβγ, is modeled using Michaelis-Men en kine ics. In [25], V2
was ini ially modeled as:
V2=k1[EP2∗][αsβγ]
[αsβγ] + K1
,
whe e k1is he Gαsac i a ion a e and K1is he dissocia ion cons an o EP2 and Gαsβγ. Howe e ,
based on [10], we simpli y he equa ion o he linea o m:
V2=k1
K1
[EP2∗][αsβγ].
3. GTP hyd olysis and deac i a ion: The deac i a ion o GαGTP
s ia GTP hyd olysis, which is no enzyme-
ca alyzed, is gi en by:
V3=k2[α∗
s],
whe e k2is he Gαshyd olysis a e.
By Assump ion 1, he concen a ion o ac i a ed EP2 ecep o s, [EP2∗], ollows a Hill equa ion. This
means [EP2∗] is de e mined by he o al concen a ion o EP2 ecep o s, EP2 o , he dissocia ion cons an ,
K2 o EP2 and PGE2, and he ini ial concen a ion o PGE2 added in he expe imen . In pa icula , [EP2∗]
a ies based on he PGE2 concen a ion:
[EP2∗] = EP2 o [PGE2]n
[PGE2]n+K2
. (4)
In [12], he Hill coe icien was assumed o be n= 1, co esponding o he classical Michaelis-Men en
eac ion. This equi alence allows he e ms ”Hill equa ion” and ”Michaelis-Men en eac ion” o be used
in e changeably when n= 1. In con as , [25] in oduced a highe Hill coe icien , n= 4, o cap u e he
h eshold beha io o EP2 o a ying ligand concen a ions, ensu ing he model e lec s he quali a i e
9

beha io obse ed in expe imen s. As shown in Figu e 4, he Hill coe icien s ongly in luences ecep o
ac i a ion, pa icula ly a low PGE2 concen a ions. Fo n= 1, ecep o ac i a ion sa u a es quickly as
PGE2 inc eases. Fo n>1, ecep o ac i a ion emains minimal a low PGE2 concen a ions, c ea ing a
h eshold e ec whe e weak signals do no igge a esponse un il a c i ical concen a ion is eached.
Figu e 4: Hill Equa ion (4) o di e en Hill coe icien s and EP2 o = 0.004µM.
Using he men ioned eac ion a es, he au onomous ODE sys em go e ning he concen a ions is:
d[αsβγ]
d =V1−V2, (5a)
d[βγ]
d =V2−V1, (5b)
d[α∗
s]
d =V2−V3, (5c)
d[αs]
d =V3−V1, (5d)
whe e as s a ed be o e:
V1: e-associa ion and dissocia ion a e o Gαsβγ,
V2: dissocia ion a e o Gαsβγ,
V3: hyd olysis a e o GTP on Gαs.
The pa ame e alues o he ODE sys em (5) a e aken om [12,25] and a e lis ed in he ollowing
able:
Pa ame e Value Uni s Desc ip ion
k15 s−1Gαsac i a ion a e
k20.07 s−1Gαshyd olysis a e
k30.7 µM−1s−1Gαsβγ associa ion a e
k418.9 ×10−3s−1Gαsβγ dissocia ion a e
K10.8 µM Dissocia ion cons an o EP2 and Gαsβγ
K20.012 µM Dissocia ion cons an o EP2 and PGE2
βγ o 0.005 µM To al Gβγ concen a ion
α o
s2.3 µM To al Gαsconcen a ion
EP2 o 0.004 µM To al EP2 concen a ion
Table 1: Pa ame e s used in he ODE sys em (5).
10
Rega ding mass conse a ion, wo ela ions can be ound:
d[αsβγ]
d +d[βγ]
d = 0 =⇒[αsβγ]+[βγ] = βγ o , (6a)
d[αsβγ]
d +d[α∗
s]
d =d[αs]
d = 0 =⇒[αsβγ]+[α∗
s]+[αs] = α o
s. (6b)
Rema k 2.1 (Exis ence and uniqueness o solu ions [21]).Conside he au onomous ODE sys em:
(u′
i( ) = i(u1, ..., un), i= 1, ..., n,
ui(0) = u0
i,i= 1, ..., n.
whe e i:D→Rna e con inuous unc ions, D⊂Rnis an open se and u(0) = (u1(0), ..., un(0)) ∈D.
I he unc ions ia e con inuously di e en iable on D, hen iis locally Lipschi z con inuous. By he
Pica d-Lindel¨o heo em, he e exis s a unique local solu ion u: [0, T)→D o he ODE sys em, whe e
T>0.
Rema k 2.2 (Non-nega i i y c i e ion).Conside he au onomous ODE sys em:
(u′
i( ) = i(u1, ..., un), i= 1, ..., n,
ui(0) = u0
i≥0, i= 1, ..., n.
The solu ions ui( ) emain non-nega i e o all ≥0 i and only i he ollowing condi ion holds o each i
and o all j=i:
i(u1, ..., 0, ..., un)≥0, whene e uj≥0, o all j=i.
Rema k 2.3 (Exis ence and uniqueness o solu ions o Sys em (5)).The unc ions on he igh -hand sides
o he equa ions a e con inuously di e en iable wi h espec o [αsβγ], [βγ], [α∗
s], [αs], as hey consis o
polynomial e ms. The e o e, by Rema k 2.1, o gi en ini ial concen a ions, Sys em (5) has a unique local
solu ion a ound he ini ial condi ions.
Rema k 2.4 (Non-nega i i y and boundedness o solu ions o Sys em (5)).I [αsβγ] = 0 and [βγ], [α∗
s], [αs]≥
0, he sys em sa is ies d[αsβγ]
d ≥0. The same holds o he o he ODEs. Thus, by he non-nega i i y
c i e ion in Rema k 2.2, all concen a ions [αsβγ], [βγ], [α∗
s], [αs] emain non-nega i e o all ≥0. Fu -
he mo e, he mass conse a ion ela ions (6) gua an ee ha he concen a ions emain bounded by βγ o
and/o α o
s.
Using he mass conse a ion ela ions (6), we can w i e he eac ion a es V1,V2,V3in e ms o
[βγ], [α∗
s].
V1=k3[αs][βγ]−k4[αsβγ] = k3[βγ]([βγ]−[α∗
s]−βγ o +α o
s)−k4(βγ o −[βγ]),
V2=k1
K1
[EP2∗][αsβγ] = k1
K1
[EP2∗](βγ o −[βγ]),
V3=k2[α∗
s].
The sys em o ou ODEs in (5) can be educed o he sys em o wo ODEs:
d[βγ]
d =V2−V1= (k4+k1
K1
[EP2∗])βγ o −k1
K1
[EP2∗] + k4+k3α o
s−k3βγ o [βγ] + k3[βγ][α∗
s]−k3[βγ]2,
d[α∗
s]
d =V2−V3=k1
K1
[EP2∗]βγ o −k1
K1
[EP2∗][βγ]−k2[α∗
s].
A e sol ing he sys em o he wo a iables [βγ] and [α∗
s], he concen a ions o [αsβγ] and [αs] a e
de e mined using he mass conse a ion ela ions (6).
11
2.2 cAMP p oduc ion
To model cAMP p oduc ion, we accoun o bo h he syn hesis by he enzyme AC and he deg ada ion by
he PDE enzyme, as men ioned in Sec ion 1. We conside wo di e en equa ions aken om [25] and [7].
Equa ion om [25]
To cap u e he posi i e in luence o AC on cAMP le els, we assume ha one dominan AC iso o m is
esponsible o cAMP p oduc ion. A ac ion o he o al AC, AC o , is ac i a ed by binding o Gα∗
s. We
model i using Hill’s equa ion wi h n= 1:
[AC∗] = AC o [α∗
s]
[α∗
s] + K4
.
Some ac ion o [AC∗] binds u he o Gβγ, o ming a ”supe -ac i a ed” AC∗-βγ complex wi h supe -
ac i a ion ac o C1:
[AC∗-βγ]=[AC∗]C1[βγ]
[βγ] + K6
.
The unbound ac ion o [AC∗] is gi en by:
[AC∗
unbound] = [AC∗]−[AC∗-βγ]=[AC∗]K6
[βγ] + K6
.
To model he cAMP deg ada ion ia PDE enzymes, we use Hill’s equa ion wi h n= 1 o desc ibe he
concen a ion o he PDE-cAMP complex:
[PDE-cAMP] = PDE o [cAMP]
[cAMP] + Kd
,
whe e we ocus on he PDE3 and PDE4 iso o ms, as hey play key oles in egula ing cAMP le els [12].
Thus, he equa ion desc ibing he cAMP concen a ion, sligh ly modi ied om [12] and as gi en in [25] is:
d[cAMP]
d =k9
AC o [α∗
s]
[α∗
s] + K4
C1[βγ] + K6
[βγ] + K6
−k11
PDE4 o [cAMP]
[cAMP] + K7
−k12
PDE3 o [cAMP]
[cAMP] + K8
. (7)
The pa ame e alues a e gi en in Table 2.
Equa ion om [7]
Unlike Equa ion (7), he supe -ac i a ion e m wi h ac o C1is no included he e. As a esul , he
p oduc ion e m o cAMP is gi en by:
A o [α∗
s]
[α∗
s] + Kas
.
The deg ada ion o cAMP, ca alyzed by he PDE enzyme, is modeled using a Michaelis-Men en eac ion
e m:
P o [cAMP]
[cAMP] + Kw
.
12
Pa ame e Value Uni s Desc ip ion
k96.713 s−1Ac i e cAMP p oduc ion a e
k11 0.72 s−1cAMP deg ada ion a e by PDE4
k12 6.85 s−1cAMP deg ada ion a e by PDE3
K40.2 µM Dissocia ion cons an o AC and Gαs
K60.09 µM Dissocia ion cons an o AC and βγ
K72.6 µM Dissocia ion cons an o PDE4 and cAMP
K80.15 µM Dissocia ion cons an o PDE3 and cAMP
AC o 0.029 µM To al AC concen a ion
PDE4 o 0.115 µM To al PDE4 concen a ion
PDE3 o 0.0025 µM To al PDE3 concen a ion
C111 βγ supe -ac i a ion ac o
Table 2: Pa ame e alues o cAMP p oduc ion in Equa ion (7), aken om [12].
Hence, he cAMP p oduc ion is desc ibed by:
d[cAMP]
d =kw
A o [α∗
s]
[α∗
s] + Kas
−dw
P o [cAMP]
[cAMP] + Kw
. (8)
He e, we only conside he main enzymes AC and PDE, which a e esponsible o cAMP p oduc ion and
deg ada ion, espec i ely. Fo simplici y, we assume ha a single dominan iso o m exis s o each enzyme.
This assump ion is made because we only ha e quali a i e expe imen al da a and no de ailed in o ma ion
abou he speci ic enzymes and iso o ms, which would be necessa y o a mo e e ined modeling app oach.
The pa ame e alues a e gi en in he ollowing able: No e ha since ela i e concen a ions a e
Pa ame e Value Uni s Desc ip ion
kw6.713 s−1Ac i e cAMP p oduc ion a e
Kas 0.2 µM Dissocia ion cons an o AC6 and Gαs
dw8.66 s−1cAMP deg ada ion a e by PDE4
Kw1.21 µM Dissocia ion cons an o PDE4 and cAMP
A o 0.0497 µM To al AC concen a ion
P o 0.039 µM To al PDE concen a ion
Table 3: Pa ame e alues o cAMP p oduc ion in Equa ion (8), aken om [7].
measu ed in he expe imen s, we assume ha he e is no basal cAMP p oduc ion a e in Equa ions (7)
and (8).
13
2.3 ODE model o cAMP p oduc ion in EP2 ecep o signaling
Combining he models o he G-p o ein ac i a ion cycle and cAMP p oduc ion, we ob ain:
d[βγ]
d = (k4+k1
K1
[EP2∗])βγ o −k1
K1
[EP2∗] + k4+k3α o
s−k3βγ o [βγ] + k3[βγ][α∗
s]−k3[βγ]2, (9a)
d[α∗
s]
d =k1
K1
[EP2∗]βγ o −K1[EP2∗][βγ]−k2[α∗
s], (9b)
d[cAMP]
d =kw
[α∗
s]A o
[α∗
s] + Kas
−dw
P o [cAMP]
[cAMP] + Kw
. (9c)
whe e [EP2∗] is gi en, as in (4), by [EP2∗] = EP2 o [PGE2]n
K3+[PGE2]n,n= 4, and he pa ame e s a e gi en in Tables
1and 3.
Rema k 2.5 (Exis ence and uniqueness, non-nega i i y and boundedness o solu ions o Sys em (9)).The
equa ions o [βγ] and [α∗
s] a e uncoupled om he one o [cAMP]. By Rema ks 2.3 and 2.4, Equa ions
(9a) and (9b) ha e a unique local solu ion ha is non-nega i e and bounded i he ini ial condi ions sa is y
0≤[βγ]0≤βγ o , 0 ≤[α∗
s]0≤α o
s.
Since he eac ion e ms in he ODE o [cAMP] a e con inuously di e en iable (Kas and Kwa e
posi i e), Rema k 2.1 gua an ees he exis ence o a unique local solu ion o Equa ion (9c). Fu he mo e,
by Rema k 2.2, he cAMP concen a ion emains non-nega i e. Howe e , he solu ion o [cAMP] is no
necessa ily bounded.
2.3.1 Simula ions
Ini ial condi ions
The ini ial condi ions o [βγ], [α∗
s] a e se o hei basal s eady-s a e concen a ions. Following [12], we
compu e he s eady-s a es o he sys em in he absence o he ligand. To de e mine hese s eady-s a es,
we se [PGE2] = 0, which using Equa ion (4), leads o [EP2∗] = 0. A he s eady-s a e, he eac ion a es
balance, so V1=V2=V3. Unde hese condi ions we ind [α∗
s]0= 0, and o [βγ] we ob ain a quad a ic
equa ion:
k3[βγ]2+ (k4+k3α o
s−k3βγ o )[βγ]−k4βγ o = 0
This equa ion gi es a single biologically ele an posi i e solu ion: [βγ]0= 5.81 ×10−5≈6×10−5.
Fu he mo e, we assume an ini ial concen a ion o ze o o cAMP, [cAMP]0= 0. Fo he nume ical sim-
ula ions, we use he odein unc ion in Py hon, see Appendix A.
Figu e 5c and 5d show ha equa ions (7) and (8) ha e di e en s eady-s a e concen a ions o high
ini ial condi ions o [PGE2]. Howe e , bo h equa ions exhibi simila quali a i e beha io , pa icula ly in
modeling he cAMP p oduc ion h eshold be ween low and high ligand concen a ions.
Gi en he lack o de ailed expe imen al da a on speci ic enzymes and iso o ms, we use Equa ion (8) o
i s educed complexi y. The simpli ying assump ions in Equa ion (8) include a single dominan iso o m, he
omission o he supe -ac i a ion e m and ewe pa ame e s. These assump ions make he model easie o
analyze bo h analy ically and nume ically, while s ill cap u ing he essen ial h eshold beha io o cAMP
p oduc ion.
14

(a) Gβγ concen a ion. (b) Gαsconcen a ion.
(c) cAMP p oduc ion using Equa ion (7). (d) cAMP p oduc ion using Equa ion (8).
Figu e 5: Solu ions o Sys em (9).
In [25], he h eshold be ween low and high PGE2 concen a ions was modeled by inc easing he Hill
coe icien o n= 4 in Equa ion (4), which gi es he concen a ion o ac i a ed ecep o s [EP2∗]. A Hill
coe icien g ea e han 1 indica es posi i e coope a i i y, ypically sugges ing ha he ecep o has mul iple
ligand-binding si es. Howe e , GPCRs gene ally ha e only a single binding si e [25]. This coope a i i y
could po en ially a ise om ecep o dime iza ion. While [25] has a mo e de ailed discussion on his h esh-
old mechanism and explo es possible mechanisms o his beha io , no de ini i e biological a gumen can be
gi en o he Hill coe icien n= 4. Fu he expe imen s a e needed o de e mine which Hill exponen bes
i s he expe imen al da a o c i ical ligand concen a ions nea he h eshold, speci ically [PGE2] ∈0.1, 1.
Thus, Sys em (9) success ully ep oduces he quali a i e expe imen al dynamics shown in Figu e 3. A
low ligand concen a ions, he EP2 ecep o emains inac i e, esul ing in no cAMP p oduc ion. Howe e ,
a highe ligand concen a ions, cAMP p oduc ion inc eases and s abilizes a a s eady-s a e le el.
15
2.4 Modeling ligand- ecep o dynamics
In [25] and he p e ious sec ion, he ligand- ecep o in e ac ion is assumed o be apid allowing he use
o a quasi-s eady-s a e app oxima ion o he ac i a ed ecep o concen a ion (Equa ion (4)). He e, we
aim o ex end he model and include he ull dynamics o he ligand- ecep o in e ac ion o e alua e he
alidi y o his assump ion.
Ou app oach is based on he model p esen ed in [9], which ma hema ically o mula es ligand- ecep o
in e ac ion-d i en cell mig a ion in he p esence o decoy ecep o s. By inco po a ing he ligand- ecep o
dynamics in o Sys em (9), we simula e he ull model and compa e he o iginal model (9), which uses
he s eady-s a e app oxima ion o [EP2∗], wi h he e o mula ed model ha explici ly includes he ligand-
ecep o dynamics.
Acco ding o [7], EP2 ecep o s do no in e nalize, so we se he in e naliza ion a e ki= 0 o ou EP2
ecep o s, see [9].
The ac i a ion o he EP2 ecep o , ep esen ed by EP2∗, occu s h ough ligand- ecep o binding wi h
he ligand PGE2. This p ocess is desc ibed by he e e sible eac ion:
EP2 + PGE2 ka
⇌
kd
EP2∗, (10)
whe e kais he associa ion o binding a e cons an , and kdis he dissocia ion o deg ada ion a e cons an .
Using he law o mass ac ion, p e iously in oduced in Sec ion 2, we de i e he ollowing sys em o ODEs
o model he dynamics o each concen a ion in his eac ion. Again, le [·] ep esen he concen a ion o
each subs ance. Then, he modi ied equa ions om [9], go e ning he ecep o -ligand dynamics a e gi en
by:
d[PGE2]
d =U1−U2=kd[EP2∗]−ka[PGE2][EP2], (11a)
d[EP2∗]
d =U2−U1=ka[PGE2][EP2] −kd[EP2∗], (11b)
d[EP2]
d =U1−U2=kd[EP2∗]−ka[PGE2][EP2]. (11c)
Fi s , we ou line he mass conse a ion ela ions o he sys em
d[EP2∗]
d +d[EP2]
d = 0 =⇒[EP2∗] + [EP2] = EP2 o , (12a)
d[PGE2]
d +d[EP2∗]
d = 0 =⇒[PGE2] + [EP2∗] = PGE2 o . (12b)
Thus, he sys em o h ee dependen equa ions can be educed o one independen equa ion, desc ibing
he dynamics o he whole Sys em (11). Howe e , we only make use o he ac i a ed and non-ac i a ed
ecep o mass conse a ion (12a).
16
Secondly, we look o he s eady-s a es o he sys em. The ecep o -ligand sys em eaches he s eady-
s a e when he ime de i a i es in Equa ions (11) become ze o, ha is, he a es U1and U2become
equal.
U2−U1= 0 ⇐⇒ ka[PGE2][EP2]−kd[EP2∗] = 0 ⇐⇒ ka[PGE2](EP2 o −[EP2∗])−kd[EP2∗]=0 ⇐⇒
⇐⇒ (ka[PGE2] + kd)[EP2∗] = kaEP2 o [PGE2] ⇐⇒ [EP2∗] = kaEP2 o [PGE2]
ka[PGE2] + kd
=EP2 o [PGE2]
[PGE2] + kd/ka
.
Hence, he s eady-s a e exp ession o he ac i a ed EP2 ecep o is gi en by:
[EP2∗] = EP2 o [PGE2]
[PGE2] + K,
whe e K=kd
ka. This exp ession o [EP2∗] ma ches Equa ion (4), wi h n= 1. To gene alize he sys em o
any Hill coe icien n, pa icula ly n= 4, o model he h eshold beha io , we eplace [PGE2] by [PGE2]n
in Equa ions (11). This esul s in he desi ed sys em o ecep o -ligand dynamics:
d[PGE2]
d =U1−U2=kd[EP2∗]−ka[PGE2]n[EP2], (13a)
d[EP2∗]
d =U2−U1=ka[PGE2]n[EP2] −kd[EP2∗], (13b)
d[EP2]
d =U1−U2=kd[EP2∗]−ka[PGE2]n[EP2]. (13c)
To es ima e he associa ion cons an ka, we use he dissocia ion a e kd= 0.0058s−1 om [9] and he
equilib ium cons an K2=K= 0.012µM om [25], and we compu e ka=kd/K2.
The inal sys em o equa ions o he simula ion, aking in o accoun ligand- ecep o dynamics, includes
he dynamics o [βγ], [α∗
s] and [cAMP]:
d[PGE2]
d =kd[EP2∗]−ka[PGE2]n[EP2], (14a)
d[EP2∗]
d =ka[PGE2]n[EP2] −kd[EP2∗], (14b)
d[EP2]
d =kd[EP2∗]−ka[PGE2]n[EP2], (14c)
d[βγ]
d = ( k1
K1
[EP2∗] + k4)βγ o −k1
K1
[EP2∗] + k4+k3α o
s−k3βγ o [βγ] + k3[βγ][α∗
s]−k3[βγ]2, (14d)
d[α∗
s]
d =k1
K1
[EP2∗]βγ o −k1
K1
[EP2∗][βγ]−k2[α∗
s], (14e)
d[cAMP]
d =kw
[α∗
s]A o
[α∗
s] + Kas
−dw
P o [cAMP]
[cAMP] + Kw
. (14 )
Rema k 2.6 (Exis ence and uniqueness, non-nega i i y and boundedness o Sys em (14)).By Rema ks 2.1
and 2.2, Equa ions (14a), (14b) and (14c) ha e unique local solu ions ha a e non-nega i e and bounded
by EP2 o and/o PGE2 o o non-nega i e ini ial condi ions. Simila ly, by Rema k (2.5), Equa ions (14d),
(14e) and (14 ) also admi unique, non-nega i e local solu ions.
17
2.4.1 Simula ions
Ini ial condi ions
The ini ial condi ions o [βγ], [α∗
s] and [cAMP] a e he same as in Sys em (9): [α∗
s]0= [cAMP]0= 0 and
[βγ]0= 6 ×10−5.
We conside ou di e en ini ial alues o he ligand concen a ion: [PGE2]0∈ {0.01, 0.1, 1, 10}. The
ini ial concen a ion o ac i a ed ecep o s is se o ze o, while he non-ac i a ed ecep o concen a ion is
de e mined using he mass conse a ion ela ion (12). Thus, [EP2∗]0= 0 and [EP2]0= EP2 o .
The plo s in Figu e 6con i m he alidi y o he apid ligand- ecep o dynamics assump ion, jus i ying
he use o he quasi-s eady-s a e app oxima ion o [EP2∗]. The EP2 ecep o s ac i a e quickly, esul ing
in no quali a i e di e ences be ween he plo s in Figu es 5and 6. Consequen ly, Sys em (9) is su icien
o desc ibe he EP2 signaling pa hway.
Mo eo e , he ligand concen a ion [PGE2] emains nea ly cons an h oughou he p ocess. This s a-
bili y a ises because he binding o PGE2 o EP2 in ol es only a small ac ion o he o al ligand, leading
o negligible changes in [PGE2] du ing ecep o in e ac ion.
While he equa ions o [PGE2] and [EP2] a e s uc u ally iden ical (see Equa ions (14a) and (14c)),
hei ini ial condi ions di e signi ican ly. Speci ically, [EP2] s a s a 0.004, while [PGE2] akes alues o
0.01, 0.1, 1, and 10. This di e ence in ini ial concen a ions signi ican ly in luences he dynamics obse ed
in Figu e 6. This highligh s how ini ial condi ions can de e mine he beha io .
We ema k he e ha he h eshold obse ed in [25] occu s due o he nonlinea i y in oduced by he
ligand concen a ion in he s eady-s a e app oxima ion o he ac i a ed ecep o s, EP2∗. Ma hema ically,
his is cha ac e ized by he dependence o [EP2∗] on [PGE2]n. Fo small ligand concen a ions, he e m
[PGE2]nbecomes negligible, p e en ing ecep o ac i a ion. Con e sely, o su icien ly la ge ligand con-
cen a ions ([PGE2] ≫1), he e m [PGE2]ndomina es, apidly sa u a ing ecep o ac i a ion and d i ing
[EP2∗] o i s s eady-s a e alue.
18
b(k)=










u(k),n
0+ ∆ · (k),n
0+∆
∆xg(k),n+1
0
u(k),n
j+ ∆ · (k),n
j
.
.
.
u(k),n
Nx−1+ ∆ · (k),n
Nx−1+∆
∆xg(k),n+1
L










.
We hen sol e he linea sys em A(k)u(k),n+1 =b(k)a each ime s ep using spsol e in Py hon.
No e ha he e ms g(k),n+1
0and g(k),n+1
Lco espond o he luxes a he bounda ies, which in Sys em
(15), depend on he solu ion o he uncoupled ligand- ecep o ODE sys em. These bounda y luxes mus
be upda ed a each ime s ep using he solu ion o he ODE sys em. Then, his ODE sys em mus be
sol ed a each ime s ep using he Backwa d Eule implici me hod:
u(k),n+1 =u(k),n+ ∆ · (k)( n+1,u(k),n+1).
The solu ions o his ODE sys em a each ime s ep a e used o upda e he lux bounda y condi ions o
he eac ion-di usion sys em.
The o iginal ODE model, designed o ma ch he expe imen al da a, accu a ely ep esen s he dynamics
o EP2 signaling. Adding in acellula di usion o he model, while keeping he same pa ame e alues and
equa ions as in he o iginal ODE model, shows ha he key signaling dynamics emain localized a he
memb ane. Figu es 8,9and 10 as well as Sys em (15) show cAMP p oduc ion occu s on he cell mem-
b ane, wi h i s concen a ion decaying in he in acellula space. This decay, howe e , is no signi ican ly
in luenced by di usion. This sugges s ha in acellula di usion is negligible, and he simple ODE model
is su icien o desc ibe he sys em.
Addi ionally, he memb ane dynamics closely esemble hose desc ibed by he o iginal ODE model.
Tha is, Sys em 15 p ese es he swi ch-like beha io o di e en ini ial condi ions. Pa icula ly, he small
p oduc ion o cAMP obse ed in Figu e 8 e lec s bounda y condi ions in luencing [βγ] and [α∗
s], p ese ing
he swi ch-like beha io unde di e en ini ial ligand condi ions. Howe e , when adding ini ial ligand con-
cen a ions o 1µM o 10µM, he p oduc ion o cAMP on he cell memb ane inc eases, see Figu es 9and 10.
Gi en ha mos concen a ions a e localized a he bounda ies o he domain, his mo i a ed he
de elopmen o a second PDE model ocused on he 1-dimensional cell memb ane o be e explo e spa ial
e ec s in ha egion.
25

Figu e 8: Le : ini ial condi ions o sys em (15). Righ : Solu ion o (15) o [PGE2]0= 0.1µM a ime
T=100.
Figu e 9: Solu ion o (15) o [PGE2]0= 1µM a ime T=100.
Figu e 10: Solu ion o (15) o [PGE2]0= 10µM a ime T=100.
26
3.2 1D eac ion-di usion model on he cell memb ane
The eac ion-di usion model in he p e ious sec ion showed ha di usion o subs a es wi hin he in acellu-
la space is negligible. This sugges s ha he ele an signaling p ocesses a e localized a he cell memb ane,
whe e he ele an subs a es a e concen a ed. Consequen ly, we in oduce ano he eac ion-di usion
model on he simpli ied 1-dimensional domain ep esen ing he cell memb ane, deno ed by Ωc= (0, L).
This app oach is jus i ied by he memb ane localiza ion o EP2 ecep o s and adenylyl cyclase (AC), he
enzyme esponsible o cAMP syn hesis [8,13].
To inco po a e spa ial di usion along he memb ane, di usion e ms a e in oduced o each subs a e in
he sys em, including ligands, ecep o s (bo h ac i a ed and non-ac i a ed o ms), and G-p o ein subuni s.
Smalle molecules, such as cAMP and ligands, a e assumed o di use as e han la ge , memb ane-bound
s uc u es like ecep o s. The di usion coe icien s a e he same as he ones used p e iously. Hence, di -
usion occu s along he cell memb ane, allowing subs a es o mo e wi hin he domain Ωcand in e ac
h ough he de ined eac ion kine ics.
All eac ions a e assumed o happen in he domain. Tha is, Reac ions (1), (2) and (3) co esponding
o he G-p o ein cycle, ligand- ecep o binding (10) and cAMP p oduc ion and deg ada ion by enzymes AC
and PDE espec i ely, a e conside ed in he domain Ωc.
Rega ding he bounda y condi ions, we apply homogeneous Neumann bounda y condi ions on he
bounda y Γc={0, L}. This choice ep esen s he assump ion o no lux o chemicals ac oss he bounda y,
meaning ha he molecules (ligand, ecep o s and G-p o ein subuni s) do no en e o lea e he de ined
1-dimensional memb ane domain, Ωc. These bounda y condi ions a e easonable because hey e lec he
es ic ion o signaling e en s and molecules o he cell memb ane, whe e he main p ocesses o EP2 sig-
naling occu .
Al e na i ely, pe iodic bounda y condi ions can be applied, o emphasize he closed-cu e s uc u e o
he cell memb ane.
Figu e 11: 1-dimensional cell-memb ane.
27
We ake he ODE Sys em (14) and add di usion o e e y subs a e. This model is ep esen ed by he
ollowing sys em o PDEs o e he domain, Ωc:
∂[PGE2]
∂ =D1
∂2[PGE2]
∂x2+U1−U2=D1
∂2[PGE2]
∂x2+kd[EP2∗]−ka[PGE2]n[EP2], (20a)
∂[EP2∗]
∂ =D3
∂2[EP2∗]
∂x2+U2−U1=D3
∂2[EP2∗]
∂x2+ka[PGE2]n[EP2] −kd[EP2∗], (20b)
∂[EP2]
∂ =D3
∂2[EP2]
∂x2+U1−U2=D3
∂2[EP2]
∂x2+kd[EP2∗]−ka[PGE2]n[EP2], (20c)
∂[βγ]
∂ =D2
∂2[βγ]
∂x2+V2−V1=D2
∂2[βγ]
∂x2+k4βγ o +k1
K1
βγ o [EP2∗]
−(k4+k3α o
s−k3βγ o )[βγ]−k1
K1
[EP2∗][βγ] + k3[βγ][α∗
s]−k3[βγ]2, (20d)
∂[α∗
s]
∂ =D2
∂2[α∗
s]
∂x2+V2−V3=D2
∂2[α∗
s]
∂x2+k1
K1
βγ o [EP2∗]−k1
K1
[EP2∗][βγ]−k2[α∗
s], (20e)
∂[cAMP]
∂ =D1
∂2[cAMP]
∂x2+kw
[α∗
s]A o
[α∗
s] + Kas
−dw
P o [cAMP]
[cAMP] + Kw
, (20 )
wi h homogeneous Neumann bounda y condi ions on he bounda y Γc={0, L}. As in he p e ious sec ion,
we in oduce k o compac ness o no a ion:
k=(1, on Γ1
c={x= 0},
0, on Γ2
c={x=L}.
Then, he bounda y condi ions on Γca e gi en by:
(−1)kD1
∂[PGE2]
∂x= 0,
(−1)kD3
∂[EP2∗]
∂x= 0,
(−1)kD3
∂[EP2]
∂x= 0,
(−1)kD2
∂[βγ]
∂x= 0,
(−1)kD2
∂[α∗
s]
∂x= 0,
(−1)kD1
∂[cAMP]
∂x= 0.
Rema k 3.3.No e ha we conside he same di usion coe icien o he G-p o eins, allowing us o apply
he mass conse a ion ela ions (6) and educe he sys em o ou equa ions o wo equa ions. We also
assume equal di usion coe icien s o he ac i a ed and non-ac i a ed ecep o s. Howe e , we do no apply
he mass conse a ion ela ion (12) in his case.
Rema k 3.4 (Exis ence, uniqueness, non-nega i i y and boundedness o solu ions o Sys em (20)).Fo he
exis ence and uniqueness o solu ions, we e e o [26]. The non-nega i i y c i e ion in Rema k 2.2 also
applies o eac ion-di usion sys ems (see [26]). Since he eac ion e ms emain unchanged, he solu ions
emain non-nega i e.
28
3.2.1 Simula ions
Ini ial Condi ions
As in he ODE model, we assume basal ini ial concen a ions o all subs a es:
[βγ](x, 0) = 6×10−5, [α∗
s](x, 0) = 0, [EP2∗](x, 0) = 0, [PGE2](x, 0) ∈ {0.01, 0.1, 1, 10},x∈(0, L).
To e lec biological e idence ha EP2 ecep o s a e clus e ed in speci ic egions along he cell memb ane
[6], we in oduce a spa ially he e ogeneous ini ial condi ion o [EP2] wi h concen a ion peaks. We model
his dis ibu ion as a sum o Gaussian unc ions cen e ed a dis inc poin s along he domain:
[EP2](x, 0) =
N
X
i=1
hiexp −(x−pi)2
2σ2,
whe e Nis he numbe o peaks, hiis he ampli ude o he i- h peak, piis i s cen e posi ion, and σ
con ols he wid h o he peaks, ep esen ing ecep o clus e ing.
To ensu e ha he esul s o his PDE model a e compa able o hose o he ODE model, we ake
he heigh o he peaks o be equal o EP2 o = 0.004. This choice o he ini ial condi ions p o ides a
he e ogeneous and biologically ealis ic ini ial dis ibu ion o [EP2]. Mo eo e , he ini ial condi ions a e
compa able wi h he alues om he ODE model.
Implemen a ion
The disc e iza ion and implemen a ion o Sys em (20) ollows he same app oach as in he p e ious sec ion.
In he in e io poin s, he solu ion is app oxima ed using Equa ion (16). Howe e , in his case, we impose
homogeneous (ze o- lux) Neumann bounda y condi ions, which gi e:
Dk
∂u(k),n+1
j
∂xj=0 ≈Dk
u(k),n+1
0−u(k),n+1
−1
∆x= 0 ⇒u(k),n+1
−1=u(k),n+1
0,
Dk
∂u(k),n+1
j
∂xj=Nx−1≈Dk
u(k),n+1
Nx−u(k),n+1
Nx−1
∆x= 0 ⇒u(k),n+1
Nx=u(k),n+1
Nx−1.
Using hese bounda y condi ions, we can ew i e he equa ions o j= 0 and j=Nx−1 as ollows:
A j= 0 : (1 + λk)u(k),n+1
0−λku(k),n+1
1=u(k),n
0+ ∆ · (k),n
0.
A j=Nx−1 : −λku(k),n+1
Nx−2+ (1 + λk)u(k),n+1
Nx−1=u(k),n
Nx−1+ ∆ · (k),n
Nx−1.
Sys em (20) is disc e ized in o a linea sys em o equa ions:
A(k)u(k),n+1 =b(k),
whe e A(k)is he ini e di e ence ma ix (19), and b(k)is he modi ied igh -hand side ec o wi h ze o- lux
condi ions, gi en by:
b(k)
j=u(k),n
j+ ∆ · (k),n
j.
29
When simula ing he sys em wi h homogeneous ini ial condi ions (including [EP2] = EP2 o = 0.004µM),
he model exhibi s a clea h eshold beha io , like in he ODE sys em. As shown in Figu e 12, o a PGE2
concen a ion o 1µM o lowe , ecep o ac i a ion does no occu , and he e is no cAMP p oduc ion.
Howe e , when [PGE2] eaches 1µM, ecep o s become ac i a ed, and cAMP p oduc ion eaches a homo-
geneous s eady s a e along he domain.
Figu e 12: Le : PGE2 = 0.1µM. Righ : PGE2 = 1µM.
To explo e spa ial e ec s, we conside non-homogeneous ini ial condi ions o he non-ac i a ed ecep-
o s, inspi ed by e idence o ecep o clus e ing on he cell memb ane [6]. The ini ial condi ion models wo
ecep o clus e s. Again, Figu es 13,14 and 15 show a h eshold be ween low and high concen a ions
o PGE2. Howe e , nume ical simula ions show ha he di usion e ms homogenize he subs a e dis i-
bu ions o e ime, elimina ing ecep o clus e s, see Figu e 14 a T=100 and 15 a T=200. Tha is, all
concen a ions end o ha e a homogeneous dis ibu ion on he domain.
As shown in Figu es 13,14, and 15, he e is a h eshold be ween low and high concen a ions o
PGE2. Howe e , nume ical simula ions show ha he di usion e ms g adually homogenize he subs a e
dis ibu ions o e ime, dissipa ing he ecep o clus e s. This beha io is shown in Figu e 14 a T= 100
and Figu e 15 a T= 200, whe e all concen a ions end owa d a uni o m dis ibu ion ac oss he domain.
Figu e 13: Le : ini ial condi ions o Sys em (20). Righ : concen a ions o Sys em (20) o [PGE2]0=
0.1µM.
30

Figu e 14: Le : ini ial condi ions. Righ : concen a ions o Sys em (20) o [PGE2]0= 1µM a T=100.
Figu e 15: Le : ini ial condi ions. Righ : concen a ions o Sys em (20) o [PGE2]0= 1µM a T=200.
In summa y, Sys em (20) ep oduces he swi ch-like beha io obse ed in he ODE sys em, cap u ing
he h eshold-dependen ac i a ion o ecep o s and cAMP p oduc ion. Howe e , he sys em is unable o
model EP2 ecep o clus e ing along he cell memb ane. In he nex sec ion, we in oduce ligand- ecep o -
based Tu ing models o explo e hei po en ial in modeling he ecep o clus e ing on he cell memb ane.
The e o e, while bo h eac ion-di usion models (15) and (20) a e ma hema ically aluable, his ex en-
sion does no p o ide signi ican imp o emen s o e he ODE model o EP2 signaling. This sugges s on
he one hand ha in acellula di usion is negligible, making he ODE model (9) su icien , and on he
o he hand ha he second model, (20), needs o be modi ied o accoun o spa ial clus e ing.
31
4. Ligand- ecep o based Tu ing models
In he p e ious sec ion, wo eac ion-di usion models we e conside ed o s udy spa ial di usion o subs a es,
in he cy oplasm and on he cell memb ane. While hese models success ully ep oduce he h eshold beha -
io obse ed in expe imen al da a, model (20), p oposed on he cell memb ane, ailed o cap u e ecep o
clus e ing, as expe imen ally obse ed by luo escence mic oscopy. This di e ence highligh s he limi a ions
o he p e ious models in explaining spa ial pa e ns on he cell memb ane.
In his sec ion, we explo e a mechanism o p oducing spa ial pa e ns, ocusing on ecep o clus e ing
h ough ligand- ecep o in e ac ions. Fi s , we in oduce he heo y o Tu ing pa e ns, which p o ides he
unde s anding o spa ial pa e n o ma ion. Nex , we jus i y why he p e ious models could no gene a e
pa e ns. Finally, we p opose a new app oach, p oposing a class o models capable o p oducing pa e ns
unde speci ic simpli ying assump ions and condi ions. These models ha e he po en ial o be applied o
EP2 signaling when addi ional expe imen al da a becomes a ailable and may also ind b oade applica ions
in o he signaling sys ems. This wo k is inspi ed by p e ious s udies o ligand- ecep o in e ac ions in
di e en biological sys ems [11,16,15,17,27,3].
32
4.1 Gene al condi ions o di usion-d i en ins abili y
Conside he sys em o di e en ial equa ions wi h ze o lux bounda y condi ions and ini ial condi ions:















u (x, )=∆u(x, ) + γ u(x, ), (x, ),x∈Ω, >0,
(x, ) = d∆ (x, ) + γgu(x, ), (x, ),x∈Ω, >0,
(n· ∇) u
!= 0, on ∂Ω, >0,
u(x, 0) = u0(x), (x, 0) = 0(x), x∈Ω.
(21)
whe e ∂Ω is a closed bounda y o he bounded and connec ed domain Ω ⊂Rnand nis he uni ou wa d
no mal o ∂Ω. The choice o ze o lux bounda y condi ions is impo an since hese condi ions imply no
ex e nal inpu . I his was no he case, he bounda y condi ions on uand would a ec di ec ly he
spa ial pa e n [18].
We use homogeneous Neumann (ze o- lux) bounda y condi ions because hey a e ma hema ically sim-
ple o handle, biologically ele an (modeling impe meable bounda ies wi h no lux ac oss he domain),
and allow us o ocus on he sel -o ganiza ion o pa e ns wi hou ex e nal inpu . Unlike ixed (Di ichle )
bounda y condi ions, which can impose spa ial pa e ns di ec ly, ze o- lux condi ions ensu e ha any ob-
se ed spa ial s uc u es a ise solely om he in e nal dynamics o he eac ion-di usion sys em, making
hem ideal o s udying Tu ing ins abili ies and eme gen pa e ns. Howe e , Di ichle bounda y condi ions
would no g ea ly al e he ollowing analysis [18,11].
Tu ing ins abili y appea s when a eac ion-di usion sys em has a s able homogeneous s eady s a e in
he absence o di usion, which loses i s s abili y in he p esence o di usion such ha spa ial pa e ns
eme ge.
4.1.1 Linea s abili y in he absence o di usion
A spa ially homogeneous s eady-s a e (u0, 0) o he Sys em (21) sa is ies:
u =γ (u, ) = 0, =γg(u, ) = 0.
We linea ize he sys em a ound (u0, 0) by in oducing he ansla ed unc ion z= (z1,z2)Twi h z1=u−u0
and z2= − 0. Then, he linea ized sys em a ound he s eady-s a e becomes:
z =γAz,
whe e
A= u
gug (u0, 0)
= u(u0, 0) (u0, 0)
gu(u0, 0)g (u0, 0),
is he Jacobian e alua ed a he poin (u0, 0). F om now on, we w i e he pa ial de i a i es e alua ed a
he s eady-s a e wi hou hei a gumen s o simplici y. The s eady-s a e o he linea ized sys em is s able,
i Re(λ)<0 o all eigen alues o A, which o a 2-dimensional sys em is ensu ed by he condi ions:
(A) = u+g <0,
de (A) = ug − gu>0.
33
Rema k 4.1.The ace and de e minan de e mine he eigen alues and ice e sa:
(A) = λ1+λ2, de (A) = λ1λ2.
4.1.2 Di usion d i en ins abili y
Now, we aim o de i e condi ions o he sys em o become uns able unde spa ial pe u ba ions in oduced
by di usion. Conside he ull eac ion-di usion sys em (21). By applying he same pe u ba ion p ocedu e,
we linea ize he eac ion-di usion sys em a ound he equilib ium poin z= (0, 0) o ob ain he ollowing
linea ized sys em:
z =γAz + D∆z=γ u(u0, 0) (u0, 0)
gu(u0, 0)g (u0, 0)z+1 0
0d∆z. (22)
Using he me hod o sepa a ion o a iables, we look o solu ions o he o m
z(x, ) = T( )X(x),
whe e T( ) is a empo al unc ion and X(x) is a spa ial unc ion. Inse ing his ansa z in o he linea ized
sys em (22) we ge
T′( )X(x) = γAT( )X(x)+DT( )∆X(x).
Di iding by T( )X(x), we ob ain:
T′( )
T( )=γAX(x) + D∆X(x)
X(x).
He e, he le -hand side is a unc ion o only, while he igh -hand side depends on xonly. The only way
hese exp essions can be equal is i bo h a e cons an . Thus, we ha e:
T′( )
T( )=γAX(x) + D∆X(x)
X(x)=˜
λ,
o some cons an ˜
λ∈R. Hence, he sepa a ion o a iables leads o wo dis inc p oblems. Fi s , he
empo al p oblem is gi en by he equa ion
T′( ) = ˜
λT( ), (23)
which desc ibes he ime e olu ion o he sys em, and has he gene al solu ion:
T( ) = T( 0)e˜
λ ,
whe e he exponen ˜
λde e mines he empo al g ow h o decay o he solu ion.
Second, he spa ial p oblem is gi en by he ellip ic eigen alue p oblem:
(γAX(x) + D∆X(x) = ˜
λX(x), x∈Ω,
(n· ∇)X(x) = 0, x∈∂Ω.
34
•Memb ane-bound ligand- ecep o complex o ma ion is slow compa ed o he binding and u no e
dynamics. Hence, di usion is no conside ed in he equa ion o he ligand- ecep o complex C.
•Fas dynamics o he ecep o -ligand complex, compa ed o he o he chemicals. Hence, we in oduce
a quasi-s eady-s a e app oxima ion o he complex concen a ion,
−δC[C] + kon[R]m[L]n−ko [C]=0 ⇐⇒ [C] = kon
δC+ko
[R]m[L]n= Γ[R]m[L]n,
whe e Γ = kon
δC+ko . Hence, he concen a ion o bound/ac i a ed ecep o s, [C], is p opo ional o [R]m[L]n,
Assump ion 2: Complex linea dependen ecep o up egula ion.
The a e o ligand- ecep o -dependen ecep o up egula ion depends linea ly on he ligand- ecep o complex
concen a ion, [C]:
µ([C]) = [C] = Γ[R]m[L]n.
No e ha i we expec a sa u a ion o he esponse o highe ligand- ecep o concen a ions, we can
conside a Hill unc ion, howe e , his leads o a smalle Tu ing space [11],
µ([C]) = H(µ([C]), K) = H(Γ[R]m[L]n,K).
Unde hese assump ions, he model o he ligand and ecep o dynamics educes o:
˙
[L] = DL∆[L] + ρS−δL[L]−n(kon[R]m[L]n−ko [C]) =
=DL∆[L] + ρL−δL[L]−nδC[C] =
=DL∆[L] + ρL−δL[L]−nδCΓ[R]m[L]n,
˙
[R] = DR∆[R] + ρR+µ([C]) −δR[R]−m(kon[R]m[L]n−ko [C]) =
=DR∆[R] + ρR+ Γ[R]m[L]n−δR[R]−mδC[C] =
=DR∆[R] + ρR+−δR[R]+( −mδC)Γ[R]m[L]n.
These equa ions model he dynamics o ligand and ecep o concen a ions, inco po a ing eedback p o-
duc ion and ecep o -ligand in e ac ions.
In [11], i is shown ha unde speci ic assump ions, he p e ious sys em con e ges o he Schnakenbe g
model. These assump ions a e as ollows:
•Recep o -independen ligand deg ada ion is negligible compa ed o ecep o -dependen deg ada ion:
δL[L]<< nδcΓ[R]m[Ln=⇒δL= 0.
•The s oichiome y o he ligand- ecep o in e ac ion in [11] is m= 2, n= 1, co esponding o one
ecep o dime binding wo monome ic ligands. Howe e , he e, we conside gene al in ege s m,n, o
de e mine he possible combina ions.
•In [11], he eedback p oduc ion linea ly depends on he complex concen a ion wi h coe icien
= (m+n)δC. Howe e , he e, we conside wo cases o in es iga e he impac o eedback
p oduc ion:
=(0, i µ([C]) = 0
(m+n)δC, i µ([C]) = [R]m[L]n
This allows o linea eedback p oduc ion o no eedback a all.
41

We ede ine he a iables and pa ame e s o simpli y no a ion. Le U=U(X,τ) be he ecep o concen a-
ion and V=V(X,τ) he ligand concen a ion. Mo eo e , we conside he ollowing escaled pa ame e s:
ρR=k1,δR=k2,k5= Γ, k3=δCΓ, k4=ρL,δL=k6= 0.
Wi h his new no a ion, he sys em becomes:
∂U
∂τ =DU∆U+k1−k2U+ (k5−mk3)UmVn,
∂V
∂τ =DV∆V+k4−nk3UmVn.
As men ioned, we dis inguish he ollowing wo cases:
•Gene alized Schnakenbe g model
We conside linea ly complex dependen ecep o up egula ion as in [11]: µ([C]) = [R]m[L]n.
Then, = (m+n)δc, which implies, k5= (n+m)k3.
Hence, he sys em esul s in a gene alized Schnakenbe g model:
∂U
∂τ =DU∆U+k1−k2U+nk3UmVn,
∂V
∂τ =DV∆V+k4−nk3UmVn.
(28)
•No eedback p oduc ion
We do no conside complex dependen ecep o p oduc ion: µ([C]) = 0.
He e, = 0, so, k5= 0. The sys em becomes:
∂U
∂τ =DU∆U+k1−k2U−mk3UmVn,
∂V
∂τ =DV∆V+k4−nk3UmVn.
(29)
4.3.1 Gene alized Schnakenbe g model
Fi s , o ha e a model independen o he uni sys em and ha e a educed numbe o pa ame e s, we
ew i e Sys em (28) using he ollowing dimensionless a iables and pa ame e s:
u=Unk3
k21
m+n−1, =Vnk3
k21
m+n−1, =DUτ
L2,x=X
L,
d=DV
DU
,a=k1
k2nk3
k21
m+n−1,b=k4
k2nk3
k21
m+n−1,γ=L2k2
DU
.
(30)
The le -hand sides o he equa ions in Sys em (28) emain:
∂U
∂τ =∂
∂τ "k2
nk31
m+n−1
u#=k2
nk31
m+n−1∂u
∂
∂
∂τ =DU
L2k2
nk31
m+n−1∂u
∂ ,
∂V
∂τ =∂
∂τ "k2
nk31
m+n−1
#=k2
nk31
m+n−1∂
∂
∂
∂τ =DU
L2k2
nk31
m+n−1∂
∂ .
42
The Laplacians a e gi en by:
DU∆U=DU
∂2U
∂X2=DU
∂
∂X∂U
∂X=DU
∂
∂X"k2
nk31
m+n−1∂u
∂x
∂x
∂X#=DU
L2k2
nk31
m+n−1∂2u
∂x2,
DV∆V=DV
∂2V
∂X2=DV
∂
∂X∂V
∂X=DV
∂
∂X"k2
nk31
m+n−1∂
∂x
∂x
∂X#=DV
L2k2
nk31
m+n−1∂2
∂x2.
The eac ion e ms a e as ollows:
(U,V) = k2
nk31
m+n−1
u,k2
nk31
m+n−1
!=k1−k2k2
nk31
m+n−1
u+nk3k2
nk3m
m+n−1k2
nk3n
m+n−1
um n,
g(U,V) = g k2
nk31
m+n−1
u,k2
k51
m+n−1
!=k4−nk3k2
nk3m
m+n−1k2
nk3n
m+n−1
um n.
Then, we ew i e he equa ions in he dimensionless o m:
DU
L2k2
nk31
m+n−1∂u
∂ =DU
L2k2
nk31
m+n−1∂2u
∂x2+k1−k2k2
nk31
m+n−1
u−nk3k2
nk3m
m+n−1k2
nk3n
m+n−1
um n,
DU
L2k2
nk31
m+n−1∂
∂ =DV
L2k2
nk31
m+n−1∂2
∂x2+k4−nk3k2
nk3m
m+n−1k2
nk3n
m+n−1
um n.
Simpli ying he e ms we ob ain:
∂u
∂ =∂2u
∂x2+L2
DUnk3
k21
m+n−1
k1−L2
DU
k2u−L2
DU
k2um n,
∂
∂ =DV
DU
∂2
∂x2+L2
DUnk3
k21
m+n−1
k4−L2
DU
k2um n.
Hence, we ge he dimensionless sys em o he ligand- ecep o dynamics:
∂u
∂ = ∆u+γ(a−u+um n),
∂
∂ =d∆ +γ(b−um n),
(31)
whe e a,b,d,γ > 0. This sys em gene alizes he Schnakenbe g model o a bi a y s oichiome ic exponen s
mand n. No e ha o m= 2, n= 1, we eco e he well-known Schnakenbe g model, wi h eac ion e ms:
(u, ) = a−u+u2 ,g(u, ) = b−u2 .
Rema k 4.3 (Exis ence, uniqueness and non-nega i i y o solu ions o Sys em (31)).Fo he exis ence and
uniqueness o solu ions, we e e o [26]. F om he gene aliza ion o Rema k 2.2 (see [26]), solu ions o
he gene alized Schnakenbe g model emain non-nega i e.
Pa e n o ma ion
The goal is o de e mine he combina ions o s oichiome ic exponen s in he ligand- ecep o binding eac-
ion (27) ha lead o pa e n o ma ion. To do his, we ollow Tu ing’s heo y as p esen ed in Sec ion 4.1,
pe o ming a linea s abili y analysis i s in he absence o di usion and hen inco po a ing di usion-d i en
43
ins abili ies.
Homogeneous s eady-s a es
Conside Sys em (31) in he absence o di usion:
u =γ (u, ) = γ(a−u+um n),
=γg(u, ) = γ(b−um n).
The homogeneous s eady-s a es (u∗, ∗) a e de e mined by sol ing:
(u, ) = 0, g(u, ) = 0.
The esul ing homogeneous s eady-s a es a e:
(u∗, ∗) = a+b,b
(a+b)m1/n!. (32)
These homogeneous s eady-s a es, a e posi i e and eal o any a,b>0, ensu ing hei biological ele ance
in he con ex o he model.
S abili y o he s eady-s a es and Tu ing analysis
To analyze he s abili y o he s eady-s a es, we compu e he Jacobian ma ix e alua ed a (u∗, ∗):
J(u∗, ∗) =  u(u, ) (u, )
gu(u, )g (u, )(u∗, ∗)
=−1 + mum−1 nnum n−1
−mum−1 n−num n−1(u∗, ∗)
,
whe e:
u(u∗, ∗) = −1 + mb
a+b=(m−1)b−a
a+b, (u∗, ∗) = n(a+b)m
nbn−1
n>0,
gu(u∗, ∗) = −mb
a+b<0, g (u∗, ∗) = −n(a+b)m
nbn−1
n<0.
To show s abili y in he absence o di usion and di usion-d i en ins abili y, we use he necessa y
condi ions o Tu ing ins abili ies de i ed in Sec ion 4.1.2.
1. u+g <0
Subs i u ing he exp essions o uand g a (u∗, ∗), we ha e:
u+g =−1 + mb
a+b−n(a+b)m
nbn−1
n<0.
Rea anging, his gi es:
mb
a+b<1 + n(a+b)m
nbn−1
no equi alen ly (m−1)b−a<n(a+b)m
n+1bn−1
n.
2. ug − gu>0
Expanding he de e minan a (u∗, ∗), we ind:
ug − gu=−−1 + mb
a+bn(a+b)m
nbn−1
n+mb
a+bn(a+b)m
nbn−1
n.
44
Simpli ying, his educes o:
ug − gu=n(a+b)m
nbn−1
n>0,
which is always sa is ied, since a,b>0.
3. d u+g >0
Subs i u ing he de i a i es a (u∗, ∗), we ob ain:
d u+g =d−1 + mb
a+b−n(a+b)m
nbn−1
n>0.
Rea anging, his gi es:
dmb
a+b>d+n(a+b)m
nbn−1
no equi alen ly d(m−1)b−da >n(a+b)m
n+1bn−1
n.
4. d= 1
5. ug <0
Expanding ug a (u∗, ∗), we ind:
ug =−−1 + mb
a+bn(a+b)m
nbn−1
n<0.
This implies:
1−mb
a+b<0⇐⇒ 1<mb
a+b⇐⇒ a<(m−1)b,
equi ing m>1, in o de o ha e a>0.
6. (d u+g )2−4d( ug − gu)>0
E alua ing a (u∗, ∗),
d−1 + mb
a+b−nb n−1
n(a+b)m
n2
−4dnb n−1
n(a+b)m
n>0.
Summa izing, o Tu ing ins abili y o occu , he ollowing condi ions mus be sa is ied:
1. u+g <0 =⇒(m−1)b−a<n(a+b)m
n+1bn−1
n
2. ug − gu>0 (always sa is ied)
3. d u+g >0 =⇒d(m−1)b−da >n(a+b)m
n+1bn−1
n
4. d= 1
5. ug <0 =⇒a<(m−1)b=⇒m>1
6. (d u+g )2−4d( ug − gu)>0 =⇒d(−1 + mb
a+b)−nb n−1
n(a+b)m
n2>4dnb n−1
n(a+b)m
n
45
To de e mine which pai s o s oichiome ic numbe s (m,n) gi e ise o Tu ing pa e ns in Sys em (26),
we s a by no ing ha he condi ion 5implies ha no Tu ing pa e ns can occu i m= 1. The e o e, we
ocus on cases whe e m>1.
Fo example, he pai (m= 2, n= 1), co esponding o he Schnakenbe g model, is known o sa is y
he necessa y condi ions o Tu ing pa e ns o ce ain pa ame e alues and o p oduce such pa e ns.
To sys ema ically de e mine which o he pai s (m,n) allow Tu ing pa e ns, we i s analyze he necessa y
condi ions in 4.3.1 o anges o pa ame e s and hen pe o m nume ical simula ions.
He e, we conside all 12 combina ions o m∈ {2, 3, 4}and n∈ {1, 2, 3, 4}. Fo each combina ion, we
check he necessa y Tu ing condi ions 4.3.1, o e a ange o alues o he pa ame e s a∈(0, 5), b∈(0, 5)
and d∈(5, 100). The esul s show ha all hese s oichiome ic combina ions sa is y he necessa y con-
di ions o Tu ing pa e ns. Fu he mo e, nume ical simula ions con i m he eme gence o Tu ing pa e ns
o each pai , see Figu es 19,20 and 21.
Now we aim o de e mine he Tu ing space o he pa ame e s aand b, while allowing d o a y [4].
To do his, we modi y he exp essions in 4.3.1 and we w i e hem in e ms o aand u∗, gi en in (32). We
jus e o mula e condi ions 1,3and 6 om he Tu ing necessa y condi ions 4.3.1. The second condi ion is
always sa is ied, while condi ions 4and 5a e de i ed di ec ly om condi ion 3.
1. u+g <0 =⇒(m−1)(u∗−a)−a<n(u∗)m
n(u∗−a)n−1
n⇐⇒ a>u∗
m(m−1) −n(u∗)m
n(u∗−a)n−1
n,
3. d u+g >0 =⇒d(m−1)u∗−dma >n(u∗)m
n+1(u∗−a)n−1
n⇐⇒ a<u∗
m(m−1) −n
d(u∗)m
n(u∗−a)n−1
n,
6. (d u+g )2−4d( ug − gu)>0 =⇒d(m−1) −m
u∗a−n(u∗)m
n(u∗−a)n−1
n2>4nd(u∗)m
n(u∗−a)n−1
n.
The las condi ion is quad a ic on aand we can ew i e i as:
c1(d)a2+c2(d)a+c3(d)>0,
whe e we ha e:
c1(d) = d2m2,
c2(d)=2dmu∗n(u∗)m
n(u∗−a)n−1
n−d(m−1),
c3(d)=(u∗)2d2(m−1)2−2dn(u∗)m
n(u∗−a)n−1
n(m+ 1) + n2(u∗)2m
n(u∗−a)2n−1
n.
Thus, since we ha e a pa abola opening upwa ds, he quad a ic cons ain gi es wo solu ions:
a<−c2(d)−pc2(d)2−4c1(d)c3(d)
2c1(d), o a>−c2(d) + pc2(d)2−4c1(d)c3(d)
2c1(d),
and simpli ying he exp essions, we ge :
a<u∗
m
(m−1) −n
d(u∗)m
n(u∗−a)n−1
n−2sn(u∗)m
n(u∗−a)n−1
n
d
, (33)
a>u∗
m
(m−1) −n
d(u∗)m
n(u∗−a)n−1
n+ 2sn(u∗)m
n(u∗−a)n−1
n
d
. (34)
46

No e ha since d>1, condi ions 4.3.1 and 4.3.1 do no con adic . Mo eo e , he cons ain in (33)
makes he condi ion 4.3.1 supe luous. Howe e , (34) con adic s wi h condi ion 4.3.1, and canno be
sa is ied. The e o e, we ha e wo bounda y cu es (lowe and uppe ) on he (b,a) plane which bound he
Tu ing space.
a>u∗
m(m−1) −n(u∗)m
n(u∗−a)n−1
n,
a<u∗
m
(m−1) −n
d(u∗)m
n(u∗−a)n−1
n−2sn(u∗)m
n(u∗−a)n−1
n
d
.
I is easy o plo he abo e exp ession when n= 1 and m∈ {2, 3, 4}. Using u∗as a pa ame e going om
0 o ∞and using b=u∗−a.
alowe >u∗
m((m−1) −(u∗)m) , (35)
auppe <u∗
m (m−1) −(u∗)m
d−2 (u∗)m
d!. (36)
Howe e , when n∈ {2, 3, 4}, a he han using he analy ic exp essions, we can nume ically e alua e he
code used o e i y he Tu ing necessa y condi ions ac oss a ange o alues. The esul ing pai s (b,a) ha
sa is y he condi ions a e plo ed, see Figu e 18.
The lowe cu e is independen o d, hus, he size o he Tu ing space can be changed by uning d.
The size o he Tu ing space is p opo ional o d, inc easing dleads o a g ea e Tu ing space. Figu es 17
and 18 show ha b>ain mos cases, as he Tu ing space mainly lies below he line a=b. Mo eo e ,
while he shape o he Tu ing space emains simila , i s size inc eases o highe pai s o s oichiome ic
numbe s, see Figu e 18.
Figu e 17: Tu ing space o pa ame e s (b,a) o he gene alized Schnakenbe g model (31) o di e en
alues o di usion a io d. The bounding cu es a e gi en by (35), (36).
47
Figu e 18: Tu ing space o pa ame e s (b,a) o he gene alized Schnakenbe g model (31) o d= 40.
Simula ions
The implemen a ion o Sys em (31) ollows he app oach desc ibed in Sec ion 3.2. Ze o- lux bounda y
condi ions a e imposed o ensu e no lux ac oss he domain bounda ies. The ini ial condi ions o he
ecep o (u) and ligand ( ) concen a ions a e se as small ampli ude andom pe u ba ions a ound he
homogeneous s eady s a e (u∗, ∗) gi en in (32). Speci ically, he ini ial condi ions o Figu es 19,20 and
21 a e gi en by:
u0= (a+b) + 0.01 · and(Nx),
0=b
um
01
n
+ 0.01 · and(Nx).
whe e Nxis he numbe o spa ial poin s, aand ba e non-nega i e cons an s and and(Nx) gene a es
uni o mly dis ibu ed andom alues.
These simula ions aim o show i he di e en s oichiome ic pa ame e s mand ncan lead o he
eme gence o Tu ing pa e ns. We use a bi a y alues ha sa is y he Tu ing necessa y condi ions 4.3.1
o a,b(gi en in he igu e cap ions below) and we se d= 40 and γ= 100. Figu es 19,20 and 21
show ha pa e ns a e o med o all combina ions o m∈ {2, 3, 4}and n∈ {1, 2, 3, 4}. Mo eo e , he
peaks in ecep o and ligand concen a ions occu al e na ely, wi h he ecep o concen a ion peaks being
signi ican ly la ge and sha pe han hose o he ligand.
48
Figu e 19: Gene alized Schnakenbe g model (31) o m= 2, n= 1, 2, 3, 4 and a= 0.125, b= 0.42
Figu e 20: Gene alized Schnakenbe g model (31) o m= 3, n= 1, 2, 3, 4 and a= 0.2, b= 1.3.
49
Figu e 21: Gene alized Schnakenbe g model (31) o m= 4, n= 1, 2, 3, 4 and a= 0.4, b= 1.
50
A. Codes
1impo numpy as np
2impo ma plo lib . pyplo as pl
3
4# Pa ame e s
5K2 = 0.012 # mu = 12 nM
6EP2_ o = 0.004 # muM = 4nM
7
8de EP2_s a (x, n):
9 e u n (x**n * EP2_ o ) / (x**n + K2)
10
11
12 x = np . linspace (0 , 1.5 , 1000)
13
14 pl . igu e ( igsize =(8 , 5) )
15 o nin [1, 2, 3, 4, 5]:
16 pl . plo (x, EP2_s a (x , n) , label = ’n␣=␣{n}’)
17
18 # Labels and i le
19 pl . i le ( "$[ ex {EP2 }^*] $␣ o ␣ di e en ␣ Hill ’s␣ coe icien s ")
20 pl . xlabel ( " $PGE2$␣($ mu␣M$)")
21 pl . ylabel ( "$[ ex {EP2 }^*] $␣($ mu␣M$)")
22 pl . legend ( i le =" Hill ’s␣ coe icien ")
23 pl . show ()
Lis ing 1: Py hon Sc ip o Figu e 4.
1impo numpy as np
2impo ma plo lib . pyplo as pl
3 om scipy . in eg a e impo odein
4
5# Pa ame e s
6k1 = 5 #s^ -1
7k2 = 0.07 #s^-1
8k3 = 0.7 # muM ^-1 s^ -1 = 0.7 e+6 M^-1 s^-1
9k4 = 18.9e -3 #s^-1
10 be a_gamma_ o = 0.005 # muM = 5nM
11 alpha_s_ o = 2.3 # muM
12 EP2_ o = 0.004 # mu = 4nM
13 K1 = 0.8 # muM
14 K2 = 0.012 # mu = 12 nM
15
16 n=4# Hill ’s coe icien
17
18 k9 = 6.713 # s^ -1
19 k11 = 0.72 # s^-1
20 k12 = 6.85 # s^-1
21 K4 = 0.2 # muM = 200 nM
22 K6 = 0.09 # muM = 90 nM
23 K7 = 2.6 # muM
24 K8 = 0.15 # muM
57

25 AC_ o = 0.029 # muM = 28 nM
26 PDE4_ o = 0.115 # muM = 115 nM
27 PDE3_ o = 0.0025 # muM = 2.5 nM
28 C1 = 11
29
30 kw = 6.713 #s^-1
31 Kas = 0.2 # muM
32 dw = 8.66 #s^-1
33 Kw = 1.21 # muM
34 A_ o = 0.0497 # muM
35 P_ o = 0.039 # muM
36
37 kd = 0.0058 # s ^-1 ( dissocia ion cons an )
38 ka = kd / K2 # Associa ion cons an
39
40 # Calcula e EP2_s a o each PGE2 alue
41 PGE2_ als = [0.01 , 0.1 , 1, 10]
42 # PGE2_ als = np. linspace (0.1 ,1 ,10)
43 EP2_s a _ als = [( PGE2 ** n * EP2_ o ) / ( PGE2 ** n + K2) o PGE2 in
PGE2_ als ]
44
45 de ODE( a s, ):
46 be a_gamma , alpha_s_s a , cAMP = a s
47 d_be a_gamma = (( k1/K1) * EP2_s a + k4) * be a_gamma_ o - (( k1/K1 ) *
EP2_s a + k4 + k3 * alpha_s_ o - k3 * be a_gamma_ o ) * be a_gamma
+ k3 * be a_gamma * alpha_s_s a - k3 * be a_gamma **2
48 d_alpha_s_s a = (k1/K1) * EP2_s a * be a_gamma_ o - (k1/K1) *
EP2_s a * be a_gamma - k2 * alpha_s_s a
49
50 # Re e ence 1
51 # d_cAMP = k9 * (( AC_ o * alpha_s_s a ) * (K6 + C1 * be a_gamma )) / ((
alpha_s_s a + K4) * (K6 + be a_gamma )) - (k11 * PDE4_ o * cAMP) / (
cAMP + K7 ) - (k12 * PDE3_ o * cAMP ) / ( cAMP + K8)
52
53 # simpli ied ( supe ac i a ion e m ou )
54 # d_cAMP = k9 * ( AC_ o * alpha_s_s a ) / ( alpha_s_s a + K4) - (k11 *
PDE4_ o * cAMP) / ( cAMP + K7) - (k12 * PDE3_ o * cAMP) / ( cAMP + K8
)
55
56 # mo e simpli ied ( linea e m )
57 # d_cAMP = k9 * ( AC_ o * alpha_s_s a ) / K4 - (k11 * PDE4_ o * cAMP ) /
(K7 ) - (k12 * PDE3_ o * cAMP ) / (K8)
58
59 # Re e ence 2
60 d_cAMP = (kw * A_ o * alpha_s_s a ) / ( alpha_s_s a + Kas) - (dw *
P_ o * cAMP ) / ( cAMP + Kw)
61
62 e u n [ d_be a_gamma , d_alpha_s_s a , d_cAMP ]
63
64 # Time ange o simula ion
65 = np . linspace (0 , 200 , 300)
66
58
67 # Ini ial condi ions
68 ic = [6e -5, 0, 0]
69
70 # Ti les and labels o he plo s
71 plo s = [
72 {"i":0," i le ": "T ajec o ies␣o ␣$[ be a gamma ]$","ylabel": "$[
be a gamma ]$␣$( mu␣ ex {M}) $"},
73 {"i":1," i le ": "T ajec o ies␣o ␣$[ alpha_s ^*] $","ylabel": "$[
alpha_s ^*]$␣$( mu ␣ ex {M }) $"},
74 {"i":2," i le ": " T ajec o ies ␣o ␣[ cAMP ]" ,"ylabel": "[ cAMP]␣$( mu ␣
ex {M })$"}
75 ]
76
77 o iin plo s :
78 pl . igu e ( igsize =(10 , 6) )
79 o PGE2 , EP2_s a in zip ( PGE2_ als , EP2_s a _ als ):
80 sol = odein (ODE , ic , )
81 pl . plo ( , sol [:, i["i"]] , label = ’PGE2␣=␣{ PGE2 :.1 }’)
82 pl . i le (i[" i le "], on size =16)
83 pl . xlabel (" Time ␣(s)", on size =14)
84 pl . ylabel (i["ylabel"], on size =14)
85 #pl . g id (T ue)
86 pl . legend ()
87 pl . show ()
88
89 de ligand_dyn ( a s, ):
90 be a_gamma , alpha_s_s a , PGE2 , EP2_s a , EP2 , cAMP = a s
91 d_be a_gamma = (( k1/K1) * EP2_s a + k4) * be a_gamma_ o - (( k1/K1 ) *
EP2_s a + k4 + k3 * alpha_s_ o - k3 * be a_gamma_ o ) * be a_gamma
+ k3 * be a_gamma * alpha_s_s a - k3 * be a_gamma **2
92
93 d_alpha_s_s a = (k1/K1) * EP2_s a * be a_gamma_ o - (k1/K1) *
EP2_s a * be a_gamma - k2 * alpha_s_s a
94
95 d_cAMP = (kw * alpha_s_s a * A_ o ) / ( alpha_s_s a + Kas ) - (dw *
P_ o * cAMP ) / ( cAMP + Kw)
96
97 # Equa ions o ligand - ecep o dynamics
98 d_PGE2 = kd * EP2_s a - ka * PGE2 **n * EP2
99 d_EP2_s a = ka * PGE2 **n * EP2 - kd * EP2_s a
100 d_EP2 = kd * EP2_s a - ka * PGE2 **n * EP2
101 e u n [ d_be a_gamma , d_alpha_s_s a , d_PGE2 , d_EP2_s a , d_EP2 , d_cAMP ]
102
103 # Time ange o he simula ion
104 = np . linspace (0 , 200 , 300)
105
106 # Ini ial PGE2 alues
107 PGE2_i = [0.01 , 0.1 , 1, 10]
108
109 # Va iables o plo and hei co esponding indices in he solu ion a ay
110 plo s = [
111 {"i":0," i le ": "T ajec o ies␣o ␣$[ be a gamma ]$","ylabel": "$[
59
be a gamma ]$␣$( mu␣M)$"},
112 {"i":1," i le ": "T ajec o ies␣o ␣$[ alpha_s ^*] $","ylabel": "$[
alpha_s ^*]$␣$( mu ␣M)$"},
113 {"i":2," i le ": " T ajec o ies ␣o ␣[ PGE2 ]" ,"ylabel": "[ PGE2]␣$( mu ␣M)
$"},
114 {"i":3," i le ": "T ajec o ies␣o ␣$[ ex { EP2 }^*]$","ylabel": "$[
ex { EP2 }^*] $␣$( mu ␣M)$"},
115 {"i":4," i le ": " T ajec o ies ␣o ␣[ EP2]" ,"ylabel": "[EP2]␣$( mu ␣M)$"
},
116 {"i":5," i le ": " T ajec o ies ␣o ␣[ cAMP ]" ,"ylabel": "[ cAMP]␣$( mu ␣M)
$"},
117 ]
118
119 # Loop o e each a iable o gene a e plo s
120 o iin plo s :
121 pl . igu e ( igsize =(10 , 6) )
122 o PGE2_ini in PGE2_i :
123 # Se he ini ial condi ions
124 ic = [6e -5, 0, PGE2_ini , 0, EP2_ o , 0]
125
126 # Sol e he sys em o ODEs
127 solu ion = odein ( ligand_dyn , ic , )
128
129 # Plo he solu ion o he cu en a iable
130 pl . plo ( , solu ion [:, i["i"]] , label = ’PGE2 ␣=␣{ PGE2_ini :.2 }’)
131
132 # Add labels , i le , and legend
133 pl . i le (i[" i le "], on size =16)
134 pl . xlabel (" Time ␣(s)", on size =14)
135 pl . ylabel (i["ylabel"], on size =14)
136 pl . legend ()
137 pl . show ()
Lis ing 2: Py hon Sc ip o Figu es 5and 6.
1impo numpy as np
2 om scipy . op imize impo sol e
3impo ma plo lib . pyplo as pl
4 om scipy . spa se impo diags
5 om scipy.spa se.linalg impo spsol e
6
7# ODE model o he ligand - ecep o dynamics o he bounda y luxes o he
PDE
8
9de backwa d_eule (PGE2 , EP2 , EP2_s a , k_d , k_a , n, d ):
10 de ligand_dyn ( a s):
11 PGE2_new , EP2_new , EP2_s a _new = a s
12 U1 = k_d * EP2_s a _new
13 U2 = k_a * ( PGE2_new **n) * EP2_new
14 e u n [
15 PGE2_new - PGE2 - d * (U1 - U2),
16 EP2_s a _new - EP2_s a - d * (U2 - U1),
17 EP2_new - EP2 - d * (U1 - U2)
60
18 ]
19
20 ini ial_condi ion = [ PGE2 , EP2 , EP2_s a ]
21 solu ion = sol e ( ligand_dyn , ini ial_condi ion )
22 e u n solu ion
23
24 # Pa ame e s
25 K2 = 0.012 # mu = 12 nM
26 k_d = 0.0058 # s^ -1 ( dissocia ion cons an )
27 k_a = k_d / K2 # s^ -1 ( associa ion cons an )
28 n=4 # Hill coe icien
29
30 d = 0.01 # Time s ep size
31 T = 200 # To al simula ion ime
32
33 num_s eps = in (T / d )
34 ime_poin s = np. linspace (0, T, num_s eps + 1)
35 PGE2_solu ions = np . ze os ( num_s eps + 1)
36 EP2_solu ions = np . ze os ( num_s eps + 1)
37 EP2_s a _solu ions = np . ze os ( num_s eps + 1)
38
39 # Ini ial condi ions
40 PGE2_0 = 10 # 0.01 ,0.1 ,1 ,10
41 PGE2_solu ions[0] = PGE2_0
42 EP2_solu ions [0] = 0.004
43 EP2_s a _solu ions [0] = 0
44
45 # Compu ing solu ions
46 o s ep in ange (0 , num_s eps ):
47 PGE2_p e = PGE2_solu ions [s ep ]
48 EP2_p e = EP2_solu ions [ s ep]
49 EP2_s a _p e = EP2_s a _solu ions [ s ep]
50 PGE2_new , EP2_new , EP2_s a _new = backwa d_eule (
51 PGE2_p e , EP2_p e , EP2_s a _p e , k_d , k_a , n, d
52 )
53 PGE2_solu ions [ s ep + 1] = PGE2_new
54 EP2_solu ions [ s ep + 1] = EP2_new
55 EP2_s a _solu ions [ s ep + 1] = EP2_s a _new
56
57 # Plo ing he solu ions
58 pl . igu e ( igsize =(10 , 6) )
59 pl . plo ( ime_poin s , EP2_solu ions , label ="EP2 ", lines yle ="-", linewid h
=2)
60 pl . plo ( ime_poin s , EP2_s a _solu ions , label =" EP2*", lines yle ="--",
linewid h =2)
61 pl . xlabel (" Time ␣(s)")
62 pl . ylabel ( " Concen a ion ␣($ mu$M)")
63 pl . i le (" T ajec o ies ␣ o ␣[ EP2 ]␣and ␣[ EP2 *]")
64 pl . legend ()
65 pl . g id ()
66 pl . show ()
67
61
68 desi ed_ ime = 10.0
69 ime_index = in (desi ed_ ime / d ) # Compu e he ime s ep index
70 EP2_s a _ alue = EP2_s a _solu ions [ ime_index ]
71
72 # Reac ion - di usion model
73
74 # Pa ame e s
75 D1 = 1
76 D2 = 0.2
77
78 k1 = 5
79 k2 = 0.07
80 k3 = 0.7
81 k4 = 18.9e -3
82 K1 = 0.8
83 be a_gamma_ o = 0.005
84 alpha_s_ o = 2.3
85 EP2_ o = 0.004
86 n=4
87
88 kw = 6.713
89 Kas = 0.2
90 dw = 8.66
91 Kw = 1.21
92 A_ o = 0.0497
93 P_ o = 0.039
94
95 # Disc e iza ion Pa ame e s
96 a=0
97 b = 70
98 T = 100 # Final ime
99 dx = 0.01
100 d = 0.01
101 Nx = in ((b - a) / dx) # Numbe o disc e iza ion poin s in space
102 N = in (T / d ) # Numbe o disc e iza ion poin s in ime
103 x = np . linspace (a , b, Nx )
104 _lin = np. a ange (0 , T + d , d )
105
106 # De ine he ma ix A ob ained om disc e iza ion
107
108 lambd_1 = D1 * (d / dx **2)
109 lambd_2 = D2 * (d / dx **2)
110
111 # Helpe unc ion o cons uc ma ices wi h ze o - lux bounda y condi ions
112 de cons uc _ma ix (Nx , lambd ):
113 uppe = -lambd * np. ones(Nx -1)
114 main = (1 + 2 * lambd ) * np. ones(Nx)
115 lowe = -lambd * np. ones(Nx -1)
116 diagonals = [lowe , main , uppe ]
117 A = diags ( diagonals , o se s =[ -1 , 0, 1] , o ma =" cs ")
118 A[0, 0] = 1 + lambd # Ze o - lux bounda y condi ion a le
119 A[Nx -1 , Nx -1] = 1 + lambd # Ze o - lux bounda y condi ion a igh
62

120 e u n A
121
122 # Cons uc ma ices
123 A_be a_gamma = cons uc _ma ix (Nx , lambd_2 )
124 A_alpha_s a = cons uc _ma ix (Nx , lambd_2 )
125 A_cAMP = cons uc _ma ix (Nx , lambd_1 )
126
127 # Sol e he sys em o equa ions using ’spsol e ’
128
129 # Ini ial condi ions
130 be a_gamma_in = np . ull (Nx , 6e -5)
131 alpha_s a _in = np. ull (Nx , 0.0)
132 cAMP_in = np. ull(Nx , 0.0)
133
134 # Lis s o s o e solu ion
135 s o e_be a_gamma = []
136 s o e_alpha_s a = []
137 s o e_cAMP = []
138
139 # Append ini ial condi ions o he lis s
140 s o e_be a_gamma . append ( be a_gamma_in . copy ())
141 s o e_alpha_s a . append ( alpha_s a _in . copy ())
142 s o e_cAMP . append ( cAMP_in . copy ())
143
144 # Time s epping
145 o ime_s ep_index , ime_s ep in enume a e ( _lin [: -1]) : # Exclude he las
ime poin
146 EP2_s a _nex = EP2_s a _solu ions [ ime_s ep_index + 1]
147
148 b_be a_gamma = be a_gamma_in + d * (
149 -k3 * be a_gamma_in * ( be a_gamma_in - alpha_s a _in -
be a_gamma_ o + alpha_s_ o )
150 + k4 * (be a_gamma_ o - be a_gamma_in))
151 # Upda e bounda ies in b ec o o be a_gamma
152 b_be a_gamma [0] += (d / dx) * (k1 / K1) * EP2_s a _nex * (
be a_gamma_ o - be a_gamma_in[0])
153 b_be a_gamma [ -1] += (d / dx) * (k1 / K1) * EP2_s a _nex * (
be a_gamma_ o - be a_gamma_in[-1])
154
155 b_alpha_s a = alpha_s a _in + d * (-k2 * alpha_s a _in )
156 # Upda e bounda ies in b ec o o alpha_s a
157 b_alpha_s a [0] += (d / dx) * (k1 / K1) * EP2_s a _nex * (
be a_gamma_ o - be a_gamma_in[0])
158 b_alpha_s a [ -1] += (d / dx) * (k1 / K1) * EP2_s a _nex * (
be a_gamma_ o - be a_gamma_in[-1])
159
160 b_cAMP = cAMP_in + d * (-dw * ( P_ o * cAMP_in ) / ( cAMP_in + Kw))
161 # Upda e bounda ies in b ec o o cAMP
162 b_cAMP [0] += kw * ( A_ o * alpha_s a _in [0]) / ( alpha_s a _in [0] + Kas )
163 b_cAMP [ -1] += kw * ( A_ o * alpha_s a _in [ -1]) / ( alpha_s a _in [ -1] +
Kas)
164
63
165 # Sol e he linea sys em Au = b o each a iable
166 be a_gamma = spsol e ( A_be a_gamma , b_be a_gamma )
167 alpha_s a = spsol e ( A_alpha_s a , b_alpha_s a )
168 cAMP = spsol e ( A_cAMP , b_cAMP )
169
170 # Upda e solu ions
171 be a_gamma_in = be a_gamma
172 alpha_s a _in = alpha_s a
173 cAMP_in = cAMP
174
175 # Upda e solu ion
176 s o e_be a_gamma . append ( be a_gamma . copy ())
177 s o e_alpha_s a . append ( alpha_s a . copy ())
178 s o e_cAMP . append ( cAMP . copy ())
179
180 # Plo solu ions o e space
181
182 ig , ax = pl . subplo s (1 , 2, igsize =(22 , 8) )
183
184 # Plo he ini ial condi ions ( ime s ep 0)
185 ax [0]. plo (x , s o e_be a_gamma [0] , ’b-’, label = ’$[ be a gamma ] _0$’,
linewid h =3)
186 ax [0]. plo (x , s o e_alpha_s a [0] , ’m-’, label = ’$[ alpha_s ^*] _0$’,
linewid h =3)
187 ax [0]. plo (x , s o e_cAMP [0] , ’y-’, label = ’[cAMP ] $_0$’, linewid h =3)
188 ax [0]. se _xlabel (’x’, on size =25)
189 ax [0]. se _ i le (" Ini ial ␣ Condi ions ", on size =20)
190
191 # Plo he inal condi ions (las ime s ep)
192 ax [1]. plo (x, s o e_be a_gamma [ -1], ’b-’, label = ’$[ be a gamma ]$(T ={T }) ’,
linewid h =3)
193 ax [1]. plo (x, s o e_alpha_s a [ -1], ’m-’, label = ’$[ alpha_s ^*] $(T ={ T}) ’,
linewid h =3)
194 ax [1]. plo (x, s o e_cAMP [ -1] , ’y-’, label = ’[ cAMP ]( T={T }) ’, linewid h =3)
195 ax [1]. se _xlabel (’x’, on size =25)
196 ax [1]. se _ i le ( " Solu ion ␣a ␣T ={T}" , on size =20)
197
198 ig . sup i le ( " Reac ion - Di usion ␣ Sys em ␣ o ␣[PGE2 ] $_0$={ PGE2_0 }␣$ mu$M",
on size =20)
199 o iin ange (2) :
200 pl . sca (ax[i])
201 pl . x icks ( on size =20 , amily =’se i ’)
202 pl . y icks ( on size =20 , amily =’se i ’)
203 ax [i ]. ick_pa ams ( axis =’bo h ’, which = ’majo ’, leng h =8)
204 ax [i ]. ick_pa ams ( axis =’bo h ’, which = ’mino ’, leng h =4)
205 ax [i]. legend ( loc =’bes ’, on size =20)
206
207 pl . show ()
Lis ing 3: Py hon Sc ip o Figu es 8,9and 10.
1impo numpy as np
2 om scipy . spa se impo diags
64
3 om scipy.spa se.linalg impo spsol e
4impo ma plo lib . pyplo as pl
5
6# Reac ion - di usion model
7# Pa ame e s
8D1 = 1
9D2 = 0.2
10 D3 = 0.08 # 0.00208 #0.2 PONER 0.05!!!
11
12 k1 = 5
13 k2 = 0.07
14 k3 = 0.7
15 k4 = 18.9e -3
16 K1 = 0.8
17 be a_gamma_ o = 0.005
18 alpha_s_ o = 2.3
19 # EP2_ o = 0.004
20 n=4
21
22 K2 = 0.012
23 kd = 0.0058
24 ka = kd / K2
25
26 kw = 6.713
27 Kas = 0.2
28 dw = 8.66
29 Kw = 1.21
30 A_ o = 0.0497
31 P_ o = 0.039
32
33 # Disc e iza ion Pa ame e s
34 a=0
35 b = 70
36 T = 100 # Final ime
37 dx = 0.01
38 d = 0.01
39 Nx = in ((b - a) / dx) # Numbe o disc e iza ion poin s in space
40 N = in (T/d ) # Numbe o disc e iza ion poin s in ime
41 x = np . linspace (a , b, Nx )
42 _lin = np. a ange (0 , T + d , d )
43
44 # De ine he A ma ix ob ained om disc e iza ion
45
46 lambd_1 = D1 * (d / dx **2)
47 lambd_2 = D2 * (d / dx **2)
48 lambd_3 = D3 * (d / dx **2)
49
50 # Helpe unc ion o cons uc ma ices wi h ze o - lux bounda y condi ions
51 de cons uc _ma ix (Nx , lambd ):
52 uppe = -lambd * np. ones(Nx -1)
53 main = (1 + 2 * lambd ) * np. ones(Nx)
54 lowe = -lambd * np. ones(Nx -1)
65
55 diagonals = [lowe , main , uppe ]
56 A = diags ( diagonals , o se s =[ -1 , 0, 1] , o ma =" cs ")
57 A[0, 0] = 1 + lambd # Ze o - lux bounda y condi ion a le
58 A[Nx -1 , Nx -1] = 1 + lambd # Ze o - lux bounda y condi ion a igh
59 e u n A
60
61 # De ine ma ices o each a iable
62 A_1 = cons uc _ma ix (Nx , lambd_1 ) # PGE2
63 A_2 = cons uc _ma ix (Nx , lambd_3 ) # EP2 *
64 A_3 = cons uc _ma ix (Nx , lambd_3 ) # EP2
65 A_4 = cons uc _ma ix (Nx , lambd_2 ) # be a_gamma
66 A_5 = cons uc _ma ix (Nx , lambd_2 ) # alpha_s_s a
67 A_6 = cons uc _ma ix (Nx , lambd_1 ) # cAMP
68
69 # Ini ial condi ion o EP2_s a
70 de ic_EP2 (x , Nx , EP2_ o , num_clus e s , sigma ):
71 EP2_in = np. ze os (Nx ) # Ini ialize wi h ze os
72 clus e _cen e s = np . linspace (x [0] , x[ -1] , num_clus e s + 2) [1: -1]
73
74 o cen e in clus e _cen e s :
75 EP2_in += EP2_ o * np. exp (-(x - cen e )**2 / (2 * sigma **2) ) # Se
peak heigh o EP2_ o
76
77 e u n EP2_in
78
79 EP2_ o = 0.004 # Peak heigh equal o EP2_ o
80 num_clus e s = 2 # Numbe o clus e s
81 sigma = 5.0 # Wid h
82
83 EP2_in = ic_EP2 (x, Nx , EP2_ o , num_clus e s , sigma)
84
85 # Sol e he sys em o equa ions using ’spsol e ’
86 # Ini ial condi ions
87 PGE2_0 = 0.1
88 PGE2_in = np. ull (Nx , PGE2_0 )
89 EP2_s a _in = np. ull (Nx , 0.0)
90 EP2_ o = 0.004
91 EP2_in = np. ull (Nx , EP2_ o )
92 be a_gamma_in = np . ull (Nx , 6e -5)
93 alpha_s a _in = np. ull (Nx , 0.0)
94 cAMP_in = np. ull(Nx , 0.0)
95
96 # Lis s o s o e solu ion
97 s o e_PGE2 = []
98 s o e_EP2_s a = []
99 s o e_EP2 = []
100 s o e_be a_gamma = []
101 s o e_alpha_s a = []
102 s o e_cAMP = []
103
104 # Append ini ial condi ions o he lis s
105 s o e_PGE2 . append ( PGE2_in )
66
79 pl . x icks ( on size =20 , amily =’se i ’)
80 pl . y icks ( on size =20 , amily =’se i ’)
81 ax [i ]. ick_pa ams ( axis =’bo h ’, which = ’majo ’, leng h =8)
82 ax [i ]. ick_pa ams ( axis =’bo h ’, which = ’mino ’, leng h =4)
83 ax [i]. legend ( loc =’bes ’, on size =20)
84
85 pl . show ()
Lis ing 8: Py hon Sc ip o Figu es 19,20 and 21.
73