# Trying to solve a set of (coupled) 1st order algebraic-differential equations [closed]

I'm trying to reproduce an anterior result on bio-physics (DOI: 10.1103/PhysRevLett.91.268101). The problem is to solve a set of coupled differential-algebraic equations.

On an attempt to do so, I've written the following code:

(*Start: Defining Constants and useful relations*)
gk = 36; gna = 120; gl = 0.3; Cm = 1; vk = -12; vna = 115; vl = 10.6; \
TauIP3 = 0.14/1000; c1 := 0.185; v1 := 6/1000; v2 :=
0.11/1000; v3 := 0.9; k3 := 0.1*1000; d1 := 0.13*1000; d2 :=
1.049*1000; d3 := 0.9434*1000; d5 := 0.08234*1000; a2 :=
0.2/(1000)^2; c0 := 2.0*1000;

rules = {AlphaN -> 0.01 (10 - V[t])/(Exp[(10 - V[t])/(10)] - 1),
BetaN -> 0.125 Exp[-V[t]/80],
AlphaM -> 0.1 (25 - V[t])/(Exp[(25 - V[t])/(10)] - 1),
BetaM -> 4 Exp[-V[t]/18], AlphaH -> 0.07 Exp[-V[t]/20],
BetaH -> 1/(Exp[(30 - V[t])/10] + 1), CaER -> (c0 - Ca[t])/c1,
mi -> IP3[t]/(IP3[t] + d1), ni -> Ca[t]/(Ca[t] + d5),
Jchannel -> c1 v1 mi^3 ni^3 q[t]^3 (Ca[t] - CaER),
Jpump -> (v3 Ca[t]^2)/(k3^2 + Ca[t]^2),
Jleak -> c1 v2 (Ca[t] - CaER),
Alphaq -> a2 d2 (IP3[t] + d1)/(IP3[t] + d3),
Betaq -> a2 Ca[t], -> Ca[t] - 196.69,
Iastro -> 2.11 Log[y] UnitStep[y - 1],
Iext -> Piecewise[{{10, 10*10^3 < t < 20*10^3}}, 0]};
(*End*)

r := 0.83600008;(*This is the parameter i want to change*)

(*Solving the problem*)
PotencialDeAcao =
NDSolve[{Cm V'[t] == -gk n[t]^4 (V[t] - vk) -
gna m[t]^3 h[t] (V[t] - vna) - gl (V[t] - vl) + Iext + Iastro,
m'[t] == AlphaM (1 - m[t]) - BetaM m[t],
n'[t] == AlphaN (1 - n[t]) - BetaN n[t],
h'[t] == AlphaH (1 - h[t]) - BetaH h[t],
IP3'[t] ==
TauIP3 (IP3Estrela - IP3[t]) + Piecewise[{{r, V[t] > 50}}],
q'[t] == Alphaq (1 - q[t]) - Betaq q[t],
Ca'[t] == -Jchannel - Jpump - Jleak, Ca[0] == 70, q[0] == 0.01,
IP3[0] == 160, V[0] == 0, m[0] == 0.05, h[0] == 0.6,
n[0] == 0.33} //. rules, {V, m, n, h, q, IP3, Ca}, {t, 0,
15*10^4}, MaxSteps -> Infinity];


On the 1st part of that code all I do is define all the constant and some useful relations of the problem.

The resulting plot should be something like this image (from the article):

> This is the results I want to get: The top graph should appear with
> r=0.2 and Iext = 10 from 100<t<140 (dont matter the t axis). The
> bottom graph should appear for r=0.8 and Iext = 10 from 100<t<110).
> PS: the time points before t=100, in this setup, doesn't matter at
> all.


EDIT:

I've found a problem with the time dependence of the variables. It was differente for the V[t] (it was ms) and for the Ca/IP3[t] it was in s. That solved the graphs started to look better. The new code is updated there. But still, the response for the case when r=0.8 and the Iext is turned for 10s still not quite right. Any ideas? It might just be a constant problem.

And another problem is the computation time. It's now takign forever to the calculation. Any performance sugestions?

-
using SolveDelayed -> True will solve until t=250 but not get rid of the messages. –  user21 Apr 1 '13 at 17:37
Dropping your "real" code here will almost always attract downvotes and reduce your opportunity to get good answers (many users will just skip your question because reading your code takes too much time). Try to isolate your problem and post a fully functional shorter toy example instead –  belisarius Apr 2 '13 at 13:33
Ok thanks for the info. Thats something I did not know :) –  Rodrigo Forti Apr 2 '13 at 13:37
Your edits sound like open-ended requests to debug your code. Such questions are off-topic and are unlikely to get much positive attention. Answerable questions are focused and will have elements of general interest. –  whuber Apr 2 '13 at 14:16
OK, in the future I'll do that somewhat different. I just thaught that posting directly my code I could be more accurate on my question. Anyway, the problem I actually have here was solved by the answer. Thats why I started a new question, which was deleted already. Anyway, sorry for any incovenience. If u want i can delete the "Edit" part. –  Rodrigo Forti Apr 2 '13 at 14:40
add comment

## closed as too localized by Jens, whuber, Oleksandr R., m_goldberg, Verbeia♦Apr 6 '13 at 7:13

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

## 1 Answer

I made a list of rules to try to clean up the (in my opinion) mess of your definitions involving the unknowns in your equation.

    rules = {AlphaN -> 0.01 (10 - V[t])/(Exp[(10 - V[t])/(10)] - 1),
BetaN -> 0.125 Exp[-V[t]/80],
AlphaM -> 0.1 (25 - V[t])/(Exp[(25 - V[t])/(10)] - 1),
BetaM -> 4 Exp[-V[t]/18],
AlphaH -> 0.07 Exp[-V[t]/20],
BetaH -> 1/(Exp[(30 - V[t])/10] + 1), CaER -> (c0 - Ca[t])/c1,
mi -> IP3[t]/(IP3[t] + d1),
ni -> Ca[t]/(Ca[t] + d5),
Jchannel -> c1 v1 mi^3 ni^3 q[t]^3 (Ca[t] - CaER),
Jpump -> (v3 Ca[t]^2)/(k3^2 + Ca[t]^2),
Jleak -> c1 v2 (Ca[t] - CaER),
Alphaq -> a2 d2 (IP3[t] + d1)/(IP3[t] + d3),
Betaq -> a2 Ca[t],
y -> Ca[t] - 0.19669,
Iastro ->
2.11 Log[y](* Piecewise[{{1,Log[y]>0}},0]}*)UnitStep[y - 1]};

PotencialDeAcao =
NDSolve[{Cm V'[t] == -gk n[t]^4 (V[t] - vk) -
gna m[t]^3 h[t] (V[t] - vna) - gl (V[t] - vl) + Iext + Iastro,
m'[t] == AlphaM (1 - m[t]) - BetaM m[t],
n'[t] == AlphaN (1 - n[t]) - BetaN n[t],
h'[t] == AlphaH (1 - h[t]) - BetaH h[t],
IP3'[t] ==
TauIP3 (IP3Estrela - IP3[t]) + r Piecewise[{{1, V[t] > 50}}],
q'[t] == Alphaq (1 - q[t]) - Betaq q[t],
Ca'[t] == -Jchannel - Jpump - Jleak, Ca[0] == 0.22, q[0] == 0.01,
IP3[0] == 0.160, V[0] == 0, m[0] == 0.05, h[0] == 0.6,
n[0] == 0.33} //. rules, {V, m, n, q, h, IP3, Ca}, {t, 0, 250},
MaxSteps -> 10^6];

Plot[{Evaluate[IP3[t] /. PotencialDeAcao],
Evaluate[Ca[t] /. PotencialDeAcao]}, {t, 0, 150}, PlotRange -> All]


-
Yes, it was a bit of a mess (I'm not a real mathematica programmer). Anyway, that figure is something I'ved managed to get earlier. And its not actually right. My question about what you've done is: why change the defintion of y? It was supposed to be a heaviside function of ln(y). But i wasn't working really well so i changed to piecewise. And why y-1? –  Rodrigo Forti Apr 1 '13 at 18:38
Nvm about the y-1 question. I get it. Thanks for the help anyway. I'll probably have to check the constant values. The initial ascent of the IP3 function should not be so fast. It should be a lot more smooth (just like the picture I linked) –  Rodrigo Forti Apr 1 '13 at 18:46
I couldn't get the picture. –  b.gatessucks Apr 1 '13 at 18:47
Try to see it here: docs.google.com/document/d/… –  Rodrigo Forti Apr 1 '13 at 18:58
Edited the post, try to see if u know how to help me please! –  Rodrigo Forti Apr 1 '13 at 22:39
add comment