(*  Fachkurs Theoretische Biophysik
  Thema Neuronale Dynamik *)

Off[General::"spell1"];
Off[General::"spell"]

(* Aufgabe 1 *)
(*
Erhöhung des Membranpotentials V => Stromfluß bei (Ca^(2+))- und (SuperPlus[K])-Kanal;
(Ca^(2+))-Reaktion instantan, (SuperPlus[K])-Reaktion verlangsamt und abhängig von Membranpotential bzgl. V3 *)

(* Aufgabe 2 
-- Fitzhugh-Modell *)
eq1:=V*(a-V)*(V-1)-W;
eq2:=b*(V-W);

a=0.25; b=1;

(* Nullklinen *)
sol1=Solve[0==eq1,W];
sol2=Solve[0==eq2,W];

Nullclines=Plot[{W/.sol1,W/.sol2},{V,-0.5,1.5},AxesLabel->{V,W},PlotRange->{{-0.5,1.5},{-0.3,0.5}}]

(* Schnittpunkt der Nullklinen = stationärer Zustand *)
sol0=Solve[{0==eq1,0==eq2},{V,W}]

sol0[[1]]

(* stationäre Zustände für beliebige a, b *)
Clear[a,b]
Solve[{0==eq1,0==eq2},{V,W}]

(* Stabilitätsanalyse *)
J[x_,y_]:={{D[eq1,x],D[eq1,y]},{D[eq2,x],D[eq2,y]}}

TableForm[J[V,W]/.sol0[[1]]]

TableForm[FullSimplify[J[V,W]/.sol0[[1]]]]

a=0.25;

TableForm[J[V,W]/.sol0[[1]]]

(* Spur ist negativ für b>-0.25 => System kann nicht instabil werden *)

Eigenvalues[J[V,W]/.sol0[[1]]]

(* Realteile negativ => System stabil *)

(* alternativ detaillierte Betrachtung mit Spur und Determinante: *)
Tr[J[V,W]/.sol0[[1]]]
Det[J[V,W]/.sol0[[1]]]

(* Determinante positiv für positive b, Spur ist negativ für b>-0.25 => System kann nicht instabil werden für b>0; keine HB, also kein periodisches Feuern möglich *)

Clear[a]
c1=ContourPlot[Det[J[V,W]/.sol0[[1]]],{a,-2,2},{b,-2,2},Contours->{0},ContourShading->False];
c2=ContourPlot[Tr[J[V,W]/.sol0[[1]]],{a,-2,2},{b,-2,2},Contours->{0},ContourShading->False,ContourStyle->RGBColor[1,0,0]];
Show[c1,c2,FrameLabel->{a,b}]

Remove["Global`*"]

(* Aufgabe 3 
-- Morris-Lecar-Modell mit I(t)=0 *)
minf[V_]:=0.5(1+Tanh[(V-V1)/V2]);
winf[V_]:=0.5(1+Tanh[(V-V3)/V4]);

V1=-1.2; V2=18; V3=2; V4=30;

Plot[minf[V],{V,-40,40}]

tau[V_]:=1/Cosh[(V-V3)/(2 V4)];

eq10:=-gCa*minf[V[t]]*(V[t]-VCa)-gK*w[t]*(V[t]-VK)-gL*(V[t]-VL);
eq11:=phi*(winf[V[t]]-w[t])/tau[V[t]];

phi=0.04;
VK=-84; VL=-60; VCa=120;
gK=8; gL=2; gCa=4.4;
Cm=20;

(* Nullklinen *)
sol01=Solve[0==eq10,w[t]];
sol02=Solve[0==eq11,w[t]];

Nullclines=Plot[{w[t]/.sol01,w[t]/.sol02},{V[t],-80,80},PlotStyle->Dashing[{0.01,0.01}]]

Show[Nullclines,PlotRange->{{-70,-50},{0,0.1}}]

(* Variation der Anfangsbedingungen für V
 => Suche nach Grenzwert für Aktionspotential *)
sol1=NDSolve[{V'[t]==eq10/Cm,V[0]==-25,w'[t]==eq11,w[0]==0},{V,w},{t,0,100}];
(*AP1=ParametricPlot[Evaluate[{V[t],w[t]}/.sol1],{t,0,100},AxesOrigin->{0,0},PlotRange->{0,0.4}]*)
AP1=ParametricPlot[Evaluate[{V[t],w[t]}/.sol1],{t,0,100},AxesOrigin->{0,0}, PlotRange->{{-70,40},{0,0.45}},AspectRatio -> 1/2]

sol2=NDSolve[{V'[t]==eq10/Cm,V[0]==-20,w'[t]==eq11,w[0]==0},{V,w},{t,0,100}];
AP2=ParametricPlot[Evaluate[{V[t],w[t]}/.sol2],{t,0,100},AxesOrigin->{0,0},PlotRange->{{-70,40},{0,0.45}},AspectRatio -> 1/2]

sol3=NDSolve[{V'[t]==eq10/Cm,V[0]==-16,w'[t]==eq11,w[0]==0},{V,w},{t,0,100}];
AP3=ParametricPlot[Evaluate[{V[t],w[t]}/.sol3],{t,0,100},AxesOrigin->{0,0},PlotRange->{{-70,40},{0,0.45}},AspectRatio -> 1/2]

sol4=NDSolve[{V'[t]==eq10/Cm,V[0]==-15,w'[t]==eq11,w[0]==0},{V,w},{t,0,100}];
AP4=ParametricPlot[Evaluate[{V[t],w[t]}/.sol4],{t,0,100},PlotRange->{{-70,40},{0,0.45}},AspectRatio -> 1/2]

sol5=NDSolve[{V'[t]==eq10/Cm,V[0]==-5,w'[t]==eq11,w[0]==0},{V,w},{t,0,100}];
AP5=ParametricPlot[Evaluate[{V[t],w[t]}/.sol5],{t,0,100}, PlotRange->{{-70,40},{0,0.45}},AspectRatio -> 1/2]

Show[Nullclines,AP1,AP2,AP3,AP4,AP5,AxesLabel->{V,w},PlotRange->{{-78,75},{-0.1,0.55}}]

(* Potentialverlauf unterhalb und oberhalb des kritischen Potentials *)
Plot[Evaluate[V[t]/.sol1],{t,0,100},PlotRange->All,PlotLabel->"Aktionspotential nicht erreicht"]
Plot[Evaluate[V[t]/.sol5],{t,0,100},PlotLabel->"Aktionspotential erreicht"]

(* (Ca^(2+))-Strom als Funktion der Zeit *)
ICa[x_]:=gCa*minf[x]*(x-VCa);

Plot[ICa[Evaluate[V[t]/.sol1]],{t,0,100},PlotRange->{{0,100},{-50,0}},AxesLabel->{t,ICa}]

Plot[ICa[Evaluate[V[t]/.sol4]],{t,0,100},AxesLabel->{t,ICa}]

Plot[ICa[Evaluate[V[t]/.sol5]],{t,0,100},AxesLabel->{t,ICa}]

(* (SuperPlus[K])-Strom als Funktion der Zeit *)
IK[x_,y_]:=gK*y*(x-VK);

Plot[IK[Evaluate[V[t]/.sol1],Evaluate[w[t]/.sol1]],{t,0,100},AxesLabel->{t,IK}]

Plot[IK[Evaluate[V[t]/.sol4],Evaluate[w[t]/.sol4]],{t,0,100},AxesLabel->{t,IK}]

Plot[IK[Evaluate[V[t]/.sol5],Evaluate[w[t]/.sol5]],{t,0,100},AxesLabel->{t,IK}]

(* Aufgabe 4 
-- Morris-Lecar-Modell mit Stimulus I(t) = (Subscript[I, 0]) *)
I0=150;

(* Oszillationen im Zeitverlauf und im Phasendiagramm *)
solOsc=NDSolve[{V'[t]==(eq10+I0)/Cm,V[0]==-25,w'[t]==eq11,w[0]==0},{V,w},{t,0,200}];
OscT=Plot[Evaluate[V[t]/.solOsc],{t,0,200},AxesLabel->{t,V}]
OscPh=ParametricPlot[Evaluate[{V[t],w[t]}/.solOsc],{t,0,200},AxesLabel->{V,w},PlotRange->{{-70,60},{0,0.6}},AspectRatio -> 1/2]

(* Aufgabe 5 
-- Abbruch der Oszillationen ("black hole") durch geeignete Parameterwahl *)
I0bh=89;

(* Oszillationen für w(0)==0 *)
solbh1=NDSolve[{V'[t]==(eq10+I0bh)/Cm,V[0]==-30,w'[t]==eq11,w[0]==0},{V,w},{t,0,500}];
bhPh1=ParametricPlot[Evaluate[{V[t],w[t]}/.solbh1],{t,0,500},AxesLabel->{V,w},PlotRange->All,PlotRange->{{-70,60},{0,0.6}},AspectRatio -> 1/2]
bhT1=Plot[Evaluate[V[t]/.solbh1],{t,0,500},AxesLabel->{t,V},PlotRange->All]

(* steady state (Fokus) für w(0)==0.12 *)
solbh2=NDSolve[{V'[t]==(eq10+I0bh)/Cm,V[0]==-30,w'[t]==eq11,w[0]==0.12},{V,w},{t,0,500}];
bhPh2=ParametricPlot[Evaluate[{V[t],w[t]}/.solbh2],{t,0,500},AxesLabel->{V,w},PlotRange->All,PlotRange->{{-70,60},{0,0.6}},AspectRatio -> 1/2]
bhT2=Plot[Evaluate[V[t]/.solbh2],{t,0,500},AxesLabel->{t,V},PlotRange->All]

(* zeitlich begrenzter Strompuls derart, daß das System vom oszillatorischen Zustand in die Nähe des stationären Zustands gebracht wird *)
Istim:=If[270<t<280,10,0]

solbh3=NDSolve[{V'[t]==(eq10+I0bh+Istim)/Cm,V[0]==-30,w'[t]==eq11,w[0]==0},{V,w},{t,0,600}];
bhPh3=ParametricPlot[Evaluate[{V[t],w[t]}/.solbh3],{t,0,600},AxesLabel->{V,w},PlotRange->All,PlotRange->{{-70,60},{0,0.6}},AspectRatio -> 1/2]
bhT3=Plot[Evaluate[V[t]/.solbh3],{t,0,600},AxesLabel->{t,V},PlotRange->All]