diffusionEquations.asv

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