-
Notifications
You must be signed in to change notification settings - Fork 0
/
diffusionEquations.asv
116 lines (99 loc) · 2.81 KB
/
diffusionEquations.asv
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
function fval = diffusionEquations(t,vec,jERM,jPM,jSOC,stencilOriginal,R,r,meanEdge,ptLst,iv)
nx = iv.nx;
btot(1:nx,1) = iv.b;%40*10^-15; %micromole/micrometer^3, total CalB concentration
betot(1:nx,1) = iv.be;
C = vec(1:nx,1);
b = vec(nx+1:2*nx,1);
Ce = vec(2*nx+1:3*nx,1);
p = vec(3*nx+1:4*nx,1);
be = vec(4*nx+1:5*nx,1);
nTrig = length(ptLst);
[jSYNPCa,jSYNPp] = inputSynapse(t,stencilOriginal,R,r,meanEdge,ptLst);
% Cytosol Ca2+
% Stencil with Neuman boundary condition
stencil = stencilOriginal;
stencil = iv.Dc*stencil;
dCdt = stencil*C+iv.kbNeg*(btot-b) - iv.kbPos*b.*C+2*r.*jERM./(R.^2 - r.^2) + 2*R.*jPM./(R.^2 - r.^2); %+ 2*R.*jSOC./(R.^2 - r.^2);
for i = 1:nTrig
dCdt(ptLst(i),1) = dCdt(ptLst(i),1) + jSYNPCa(i);
end
% CalBindin
% Stencil with Neuman boundary condition
stencil = stencilOriginal;
stencil = iv.Db*stencil;
dbdt = stencil*b+iv.kbNeg*(btot-b)-iv.kbPos*b.*C;
% ER Ca2+
% Stencil with Neuman boundary condition
stencil = stencilOriginal;
stencil = iv.Dc*stencil;
dCedt = stencil*Ce - 2*jERM./r + 2*jSOC./r + 2000*10^-4*(betot-be)-10^-5*be.*Ce;
% ER CalReticulin
% Stencil with Neuman boundary condition
stencil = stencilOriginal;
stencil = iv.Db*stencil;
dbedt = stencil*be + 2000*10^-4*(betot-be)-10^-5*be.*Ce;
% IP3
pr(1:nx,1) = iv.pr;%40*10^-18; %micromol/micrometer^3, Basal IP3 concentration
% Stencil with Neuman boundary condition
stencil = stencilOriginal;
stencil = iv.Dp*stencil;
dpdt = stencil*p-iv.kp*(p-pr);
for i = 1:nTrig
dpdt(ptLst(i),1) = dpdt(ptLst(i),1) + jSYNPp(i);
end
fval = [dCdt; dbdt; dCedt; dpdt; dbedt];
end
function [jSYNPCa,jSYNPp] = inputSynapse(t,stencilOriginal,R,r,meanEdge,ptLst)
nTrig = length(ptLst);
tStartCa = 1500;
tEndCa = 1510;
tStartp = 1500;
tEndp = 3500;
jINCa = 2.5*10^-6;
slopeCa = -jINCa/10;
jINp = 0*5*10^-6;
slopep = -jINp/2000;
jSYNPCa(1:nTrig) = 0;
jSYNPp(1:nTrig) = 0;
for i = 1:nTrig
pt = ptLst(i);
flag = isEnd(pt,stencilOriginal);
if t>tStartCa && t<tEndCa
jINCay = slopeCa*(t-tStartCa)+jINCa;
if flag
jSYNPCa(i) = 2*jINCay/meanEdge;
else
jSYNPCa(i) = 2*R(pt)*jINCay/(R(pt)^2-r(pt)^2);
end
jSYNPCa
elseif t>25000 && t<25010
jINCay = slopeCa*(t-25000)+jINCa;
if flag
jSYNPCa(i) = 2*jINCay/meanEdge;
else
jSYNPCa(i) = 2*R(pt)*jINCay/(R(pt)^2-r(pt)^2);
end
jSYNPCa
else
jSYNPCa(i) = 0;
end
if t>tStartp && t<tEndp
jINpy = slopep*(t-tStartp)+jINp;
if flag
jSYNPp(i) = 2*jINpy/meanEdge;
else
jSYNPp(i) = 2*R(pt)*jINpy/(R(pt)^2-r(pt)^2);
end
else
jSYNPp(i) = 0;
end
end
end
function flag = isEnd(pt,stencilOriginal)
% Determines if selected point is end point
if nnz(stencilOriginal(pt,:)) < 3
flag = 1;
else
flag = 0;
end
end