-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsolve_py.f95
292 lines (267 loc) · 9.54 KB
/
solve_py.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
! -*- f95 -*-
! (c) 2016 - Ilya Prokin - isprokin@gmail.com - https://sites.google.com/site/ilyaprokin
! INRIA Rhone-Alpes
! STDP model : solver of ODE system with lsoda and interface with python
subroutine lsoda_calc(y, iStatus, y0, lv, t, n)
use comp_part
implicit none
external :: FEX
integer, parameter :: NEQ=32
integer, parameter :: JT = 2
integer, parameter :: LRW = 22 + NEQ * max(16, NEQ + 9)
integer, parameter :: LIW = 20 + NEQ
real*8 :: RWORK(LRW)
integer :: IWORK(LIW), JDUM
integer :: IOPT=0, ISTATE, ITASK, iCRIT
integer, parameter :: ITOL = 1
integer, intent(in) :: n
real*8, intent(in) :: lv(165)
real*8, intent(in) :: y0(NEQ), t(n)
real*8, intent(out) :: y(NEQ,n)
integer, intent(out) :: iStatus
integer, parameter :: ncrit = 7
real*8 :: ynew(NEQ), tcrit(ncrit), tcritlast
real*8 :: t_in
integer :: i
ITASK = 1
ISTATE = 1
iCRIT=1
call parameters_get(pars, lv)
call tables_make
call stims_first_tcrit_set(pars, stims_tabs, tcrit)
tcritlast = tcrit(ncrit)+pars%stimulation%num_stim/pars%stimulation%Freq
! setting optional integration parameters
IOPT=1
RWORK(5:10)=0.0
IWORK(5:10)=0
RWORK(6)=pars%integration%HMAX
IWORK(6)=int(pars%integration%MXSTEP) ! number of interally defined steps for each integration step
ynew=y0
y(:,1)=y0
t_in=t(1)
do i=2,n
call tcrit_per_stim_set_helper(t(i), tcritlast, 1d0/pars%stimulation%Freq, ITASK, RWORK, LRW, tcrit, ncrit, iCRIT)
call DLSODA(FEX,NEQ,ynew,t_in,t(i),ITOL, pars%integration%RTOL, pars%integration%ATOL, ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JDUM,JT)
if (ISTATE<0) exit
y(:,i)=ynew
end do
iStatus = ISTATE
call tables_clean
end subroutine lsoda_calc
subroutine FEX(NEQ, T, Y, YDOT)
use comp_part
implicit none
integer, intent(in) :: NEQ
real*8, intent(in) :: T, Y(NEQ)
real*8, intent(out) :: YDOT(NEQ)
call RHS(NEQ, T, Y, YDOT)
end subroutine FEX
! --- helpers to be called from python ---
! python interface to functions from comp_part
module extras
use comp_part
implicit none
contains
subroutine CaM_conc_func(Ca_cyt, n, lv, CaM)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: Ca_cyt(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: CaM(n)
integer :: i
!f2py depend(n) Ca_cyt :: n=len(Ca_cyt)
!f2py depend(n) CaM
call parameters_get(pars, lv)
do i=1,n
CaM(i) = CaM_conc(Ca_cyt(i), pars)
end do
end subroutine CaM_conc_func
subroutine CaMKIIpho(y, n, CaMKIIpho_answer)
integer, intent(in) :: n
real*8, intent(in) :: y(n,13)
real*8, intent(out) :: CaMKIIpho_answer(n)
integer :: i
!f2py depend(n) y :: n=len(y)
!f2py depend(n) CaMKIIpho_answer
forall(i=1:n)
CaMKIIpho_answer(i)=CaMKIIpho_func(y(i,:)) ! CaMKIIpho_func(y(i,:)) - Probably wrong! Not wrong for this definition
end forall
end subroutine CaMKIIpho
subroutine CaMKIIpho_4F(y, n, CaMKIIpho_answer)
integer, intent(in) :: n
real*8, intent(in) :: y(13, n)
real*8, intent(out) :: CaMKIIpho_answer(n)
integer :: i
!f2py depend(n) y :: n=len(y)
!f2py depend(n) CaMKIIpho_answer
forall(i=1:n)
CaMKIIpho_answer(i)=CaMKIIpho_func(y(:,i))
end forall
end subroutine CaMKIIpho_4F
subroutine ical_caL13(v, h, m, calo, cali, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: v(n), h(n), m(n), calo, cali(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
call caL13_tables_make(-100d0,100d0,401,pars, caL13_tabs)
do i=1,n
answer(i) = ical_caL13_func(v(i), h(i), m(i), calo, cali(i), pars)
end do
call caL13_tables_clean(caL13_tabs)
end subroutine ical_caL13
subroutine i_TRPV1(AEA, V, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: AEA(n), V(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
do i=1,n
answer(i) = pars%TRPV1%gTRPV1 * V(i) * g_TRPV1_func(AEA(i), V(i), pars)
end do
end subroutine i_TRPV1
subroutine j_ca_TRPV1(AEA, Ca_cyt, V, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: AEA(n), Ca_cyt(n), V(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
real*8 :: G_TRPV1
integer :: i
call parameters_get(pars, lv)
do i=1,n
G_TRPV1 = g_TRPV1_func(AEA(i), V(i), pars)
answer(i) = -pars%I_to_Ca_flux%TRPV1 * pars%TRPV1%gTRPV1 * V(i) * G_TRPV1
end do
end subroutine j_ca_TRPV1
subroutine i_AMPA(V, o_AMPA, gAMPAmax, n, answer)
integer, intent(in) :: n
real*8, intent(in) :: V(n), o_AMPA(n), gAMPAmax
real*8, intent(out) :: answer(n)
integer :: i
do i=1,n
answer(i) = i_AMPA_func(V(i), o_AMPA(i), gAMPAmax)
end do
end subroutine i_AMPA
subroutine i_NMDA(V, o_NMDA, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: V(n), o_NMDA(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
call NMDA_tables_make(-100d0,100d0,401,pars, NMDA_tabs)
do i=1,n
answer(i) = pars%NMDA%gNMDA * V(i) * g_NMDA_func(V(i), o_NMDA(i))
end do
call NMDA_tables_clean(NMDA_tabs)
end subroutine i_NMDA
subroutine j_ca_NMDA(V, Ca_cyt, o_NMDA, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: V(n), Ca_cyt(n), o_NMDA(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
real*8 :: G_NMDA
integer :: i
call parameters_get(pars, lv)
call NMDA_tables_make(-100d0,100d0,401,pars, NMDA_tabs)
do i=1,n
G_NMDA=g_NMDA_func(V(i), o_NMDA(i))
answer(i) = -pars%I_to_Ca_flux%NMDA * pars%NMDA%gNMDA * V(i) * G_NMDA
end do
call NMDA_tables_clean(NMDA_tabs)
end subroutine j_ca_NMDA
subroutine JIP3R_CICR(IP3, Ca_cyt, Ca_ER, h, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: IP3(n), Ca_cyt(n), Ca_ER(n), h(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
do i=1,n
answer(i) = JIP3R_CICR_func(IP3(i), Ca_cyt(i), Ca_ER(i), h(i), pars)
end do
end subroutine JIP3R_CICR
subroutine Jserca_CICR(Ca_cyt, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: Ca_cyt(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
do i=1,n
answer(i) = Jserca_CICR_func(Ca_cyt(i), pars)
end do
end subroutine Jserca_CICR
subroutine Jleak_CICR(Ca_cyt, Ca_ER, n, lv, answer)
integer, intent(in) :: n
real*8, intent(in) :: Ca_cyt(n), Ca_ER(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: answer(n)
integer :: i
call parameters_get(pars, lv)
do i=1,n
answer(i) = Jleak_CICR_func(Ca_cyt(i), Ca_ER(i), pars)
end do
end subroutine Jleak_CICR
end module extras
!!! steady state solutions
module steady_states
use comp_part
implicit none
contains
! Ca channels
subroutine h0_m0_caL13(v, n, lv, h, m)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: v(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: h(n), m(n)
integer :: i
call parameters_get(pars, lv)
call caL13_tables_make(-100d0,100d0,401,pars, caL13_tabs)
do i=1,n
h(i) = caL13_hinf(v(i))
m(i) = caL13_minf(v(i))
end do
call caL13_tables_clean(caL13_tabs)
end subroutine h0_m0_caL13
! NMDA
subroutine o0_NMDA(Glu, n, lv, o_NMDA)
! o_NMDA - probability of the open state
implicit none
integer, intent(in) :: n
real*8, intent(in) :: Glu(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: o_NMDA(n)
call parameters_get(pars, lv)
call NMDA_tables_make(-100d0,100d0,401,pars, NMDA_tabs)
o_NMDA = pars%NMDA%Alpha*Glu/( pars%NMDA%Alpha*Glu + pars%NMDA%Beta )
call NMDA_tables_clean(NMDA_tabs)
end subroutine o0_NMDA
! CICR
subroutine h0_CICR(Ca_cyt, IP3, n, lv, h_CICR)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: Ca_cyt(n), IP3(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: h_CICR(n)
real*8 :: H1(n)
call parameters_get(pars, lv)
H1 = (pars%CICR%a2 * pars%CICR%d2 * (IP3+pars%CICR%d1)/(IP3+pars%CICR%d3))
h_CICR = H1/(H1+pars%CICR%a2*Ca_cyt)
end subroutine h0_CICR
subroutine o0_d0_AMPA(Glu, n, lv, o_AMPA0,d_AMPA0)
implicit none
integer, intent(in) :: n
real*8, intent(in) :: Glu(n)
real*8, intent(in) :: lv(165)
real*8, intent(out) :: o_AMPA0(n), d_AMPA0(n)
integer :: i
call parameters_get(pars, lv)
do i=1,n
call o_d_AMPA_IC_setup(Glu(i), pars, o_AMPA0(i), d_AMPA0(i))
end do
end subroutine o0_d0_AMPA
end module steady_states