Not enough input arguments
2 views (last 30 days)
Show older comments
h=1;
tspan=0:h:500;
yo=[6*10^4 10^6 0 10^4 0 4.07 3.33 25.75 1.45 0 0 36.2 60 0 781 105 100 1.25 13.5 20 7.4 82.5];
Nstep=size(tspan,2);
xNOL=0.75:0.05:1.5;
N=size(yo,2);
n=numel(xNOL);
t=zeros(Nstep,n);
MAPdata=[115 120 117 112 112.5 113 113 113 113 113 113 113 113];
tsimulation=0:4:48;
y=zeros(Nstep,N,n);
for k=1:n
MAP=zeros(Nstep,n);
if (tspan<=max(tsimulation))
MAP(:,k)=interp1(tsimulation,MAPdata,tspan,"linear");
else
MAP(:,k)=interp1(tsimulation,MAPdata,tspan,"linear","extrap");
end
[t(:,k),y(:,:,k),MAP,Nstep,xNOL]=ode45(@(t,y,MAP,Nstep,xNOL)ActualmodelwithPinfinity(t,y,MAP(:,k),Nstep,xNOL(k)),tspan,yo);
end
function dydt=ActualmodelwithPinfinity(t,y,MAP,Nstep,xNOL)
dydt=zeros(22,1);
% pathogen load parameters start
kpg=0.008;
kd=10^-4;
um=1.1*10^2;
kpm=2*10^-2;
sm=0.48*2*5;
a=2*10^-3;
kcNA=1000;
b=1.1*10^-3;
kcMA=100;
% pathogen load parameters start
% resting neutrophils and active neutrophils parameters start
knr=2.04;
Ns=2.5*10^6;
kNRd=1.224;
kNAd=0.364;
kNP=0.047;
xNP=100.23;
hNP=2;
kNTNF=10.92;
xNTNF=400;
kN6=10.92;
xN6=393;
kN8=10.64;
xN8=554;
xN10=48;
% resting neutrophils and active neutrophils parameters end
% resting monocytes and active monocytes parameters start
kMR=1;
Minfinty=3*10^4;
kMRd=2/3;
kMP=0.052;
xMP=4.93;
kMTNF=14;
xMTNF=723;
kM6=15.7;
xM6=783.5;
hMP=1;
kMAd=0.52;
xM10=81;
% resting monocytes and active monocytes parameters end
% TNF parameters start
kTNFN=377;
kTNFM=106;
hTNFNM=1.12;
xTNFN=100*6*10^3;
xTNFM=30*3*10^2;
xGTNF=(60/2);
xLTNF=(95*10^-2*2);
xTNF6=250;
xTNF10=235;
kTNF=0.12;
wTNF=4.07;
% TNF parameters end
% IL6 parameters start
k6N=333;
k6M=189;
h6NM=3.2;
x6N=10^5;
x6M=10^3;
k6TNF=53.89;
x6TNF=350;
x66=546;
x610=40;
xL6=(115*10^-2*2);
k6=0.32;
wIL6=3.33;
k6NO=3;
x6N0=0.08;
xG6=(55/2);
% IL6 parameters end
% IL8 parameters start
k8N=17.3*10^2;
x8N=2*10^5;
h8N=1;
k8TNF=17.46;
x8TNF=300;
x810=30;
k8=0.7;
wIL8=25.75;
% IL8 parameters end
% IL10 parameters start
k10N=0.012*5;
k10M=0.022;
x10N=10^5;
x10M=10^3;
k10TNF=3.8*10^3;
k106=4.3*10^3;
x10TNF=200;
x106=400;
k10=0.293;
wIL10=1.45;
% IL10 parameters end
% Damage parameters start
f=0.07;
T=2510;
kA=2.93;
kdp=15.51;
xdp=142.4;
ud=0.33;
% Damage parameters end
% Nitric Oxide parameters start
nNOTNF=195;
kNONM=1.2*10^-6*1.5;
kNOD=6.3*10^-5*2.5;
nNO10=150;
uno=0.025;
% Nitric Oxide parameters end
% Temperature parameters start
kTTNF=1.2;
kT6=1.6;
kT10=0.4;
kT=0.7;
nTTNF=130;
nT6=140;
nT10=40;
hT=1;
Tb=36.2;
% Temperature parameters end
% Heart rate parameters start
kH=0.46;
HM=190;
BPb=118.8;
Hb=60;
hHP=4;
nHPhBP=14.6;
nHPlBP=22.11;
nHT=0.93;
tauH=0.79;
% Heart rate parameters end
% endotoxin parameters start
kpe=112.3;
npe=100.3;
kde=0.2;
hET=2;
% endotoxin parameters end
% pain perception parameters start
kPTe=2*10^-3;
kPT=0.5;
PTb=781;
% pain perception parameters end
% glucose and glycogen parameters start
sg=109.8;
sIN=0.077;
kLG=38.23;
xING=10.5;
xTNFG=107;
xIL6G=133;
kglucgly=1.74;
kglyglu=1.827;
% glucose and glycogen parameters end
% Lactate parameters start
kGL=0.0046;
xETG=100;
uL=0.0913;
Lact0=1.25;
% Lactate parameters end
% Insulin parameters start
alpha=50;
Go=105;
xINTNF=407;
xIN6=333;
kIN=0.0417;
kING=4.9;
wINSULIN=13.5;
% Insulin parameters end
% Respiratory rate parameters start
RRb=20;
kRRlact=100;
nRRlact=1;
nTempRR=0.24;
nPaO2RR=70;
tauRR=0.1;
% Respiratory rate parameters end
% PH parameters start
kPHlact=0.01;
nPaO2Ph=20;
kpH=0.5;
PHbaseline=7.4;
Pao2b=82.5;
% PH parameters end
% PaO2 parameters start
kinfPao2=0.41;
nTempPao2=0.4*1.2;
nPHPao2=7.4;
ko2Pmetab=0.03;
% PaO2 parameters end
dydt(1)=kpg*y(1)-kd*y(1)-((sm*y(1))/(um+kpm*y(1)))-a*HU2(y(1),kcNA,1)*y(3)-b*HU2(y(1),kcMA,1)*y(5);
dydt(2)=knr*y(2)*(1-(y(2)/Ns))-kNP*HU2(y(1),xNP,hNP)*y(2)*(1+kNTNF*HU2(y(6),xNTNF,2))*(1+kN6*HU2(y(7),xN6,2))*(1+kN8*HU2(y(8),xN8,3))*HYX(y(9),xN10,2)-kNRd*y(2);
dydt(3)=kNP*HU2(y(1),xNP,hNP)*y(2)*(1+kNTNF*HU2(y(6),xNTNF,2))*(1+kN6*HU2(y(7),xN6,2))*(1+kN8*HU2(y(8),xN8,3))*HYX(y(9),xN10,2)-kNAd*y(3);
dydt(4)=kMR*(1-(y(4)/Minfinty))*y(4)-kMP*HU2(y(1),xMP,hMP)*y(4)*(1+kMTNF*HU2(y(6),xMTNF,2))*(1+kM6*HU2(y(7),xM6,2))*HYX(y(9),xM10,2)-kMRd*y(4);
dydt(5)=kMP*HU2(y(1),xMP,hMP)*y(4)*(1+kMTNF*HU2(y(6),xMTNF,2))*(1+kM6*HU2(y(7),xM6,2))*HYX(y(9),xM10,2)-kMAd*y(5);
dydt(6)=(kTNFN*HU2(y(3),xTNFN,hTNFNM)+kTNFM*HU2(y(5),xTNFM,hTNFNM))*HU2(y(16),xGTNF,2)*HYX(y(18),xLTNF,1)*HYX(y(7),xTNF6,2)*HYX(y(9),xTNF10,2)-kTNF*(y(6)-wTNF);
dydt(7)=(k6N*HU2(y(3),x6N,h6NM)+k6M*HU2(y(5),x6M,h6NM))*(k6TNF*HU2(y(6),x6TNF,2)+k6NO*HU2(y(11),x6N0,2))*HYX(y(18),xL6,1)*HU2(y(16),xG6,2)*HYX(y(7),x66,1)*HYX(y(9),x610,2)-k6*(y(7)-wIL6);
dydt(8)=(k8N*HU2(y(3),x8N,h8N)+k8TNF*HU2(y(6)-wTNF,x8TNF,3))*HYX(y(9),x810,2)-k8*(y(8)-wIL8);
dydt(9)=(k10N*HU2(y(3),x10N,4)+k10M*HU2(y(5),x10M,4))*(k10TNF*HU2(y(6),x10TNF,2)+k106*HU2(y(7),x106,4))-k10*(y(9)-wIL10);
z=f*(y(3)+y(5))-T;
if (z>0)
u=z;
else
u=0;
end
dydt(10)=(u/(1+kA*y(9)))+kdp*hV(y(1),xdp)-ud*y(10);
dydt(11)=(kNONM*(y(3)+y(5))+kNOD*y(10))*HU2(y(6)-wTNF,nNOTNF,2)*HYX(y(9)-wIL10,nNO10,1)-uno*y(11);
dydt(12)=kTTNF*HU2((y(6)-wTNF),nTTNF,hT)+kT6*HU2((y(7)-wIL6),nT6,hT)-kT10*(1-HYX((y(9)-wIL10),nT10,hT))-kT*(y(12)-Tb);
fBP=zeros(Nstep,1);
for j=1:Nstep
if(MAP(:,j)<=100)
fBP(j)=HU2(BPb-MAP(:,j),nHPlBP,hHP);
else
fBP(j)=HYX(BPb-MAP(:,j),nHPhBP,hHP);
end
dydt(13)=(Hb-y(13)+(kH*(HM-Hb)*HU2((y(12)-Tb),nHT,2)*fBP(j)))/tauH;
end
dydt(14)=kpe*HU2(y(1),npe,hET)-kde*y(14);
dydt(15)=-kPTe*y(14)*y(15)-kPT*(y(15)-PTb);
dydt(16)=(sg+kLG*(y(18)-Lact0)+(kglyglu*y(17)*(1+HU2(y(6)-wTNF,xTNFG,2))*(1+HU2(y(7)-wIL6,xIL6G,2))))*HYX(y(19)-wINSULIN,xING,1)-kglucgly*y(16)-sIN*y(16)*y(19);
dydt(17)=kglucgly*y(16)-kglyglu*y(17)*(1+HU2(y(6)-wTNF,xTNFG,2))*(1+HU2(y(7)-wIL6,xIL6G,2))*HYX(y(19)-wINSULIN,xING,1);
dydt(18)=kGL*y(16)*(1-HYX(y(11),xNOL,2))*(1+HU2(y(14),xETG,1))-uL*(y(18)-Lact0);
dydt(19)=kING*(1+HU2(y(16)-Go,alpha,2))*(1-HYX(y(6)-wTNF,xINTNF,1))*(1-HYX(y(7)-wIL6,xIN6,1))-kIN*(y(19)-wINSULIN);
dydt(20)=(RRb-y(20)+(kRRlact*(1+HU2(y(18),nRRlact,1))*HU2(y(12)-Tb,nTempRR,2)*HYX(y(22),nPaO2RR,1)))/tauRR;
dydt(21)=-kPHlact*y(18)*y(21)*(1-HYX(Pao2b-y(22),nPHPao2,1))+kpH*(PHbaseline-y(21));
dydt(22)=-kinfPao2*HU2(y(12)-Tb,nTempPao2,2)*HYX(y(21),nPaO2Ph,1)-ko2Pmetab*(y(22)-Pao2b);
end
function qw=HU2(X,n,h)
qw=(X^h)/(X^h+n^h);
end
function qz=HYX(X,n,h)
qz=(n^h)/(X^h+n^h);
end
function h=hV(Z,x)
h=(Z^6)/(Z^6+x^6);
end
0 Comments
Answers (1)
Stephen23
on 6 Dec 2024
Edited: Stephen23
on 6 Dec 2024
Your function parameterisation is incorrect. This anonymous function:
@(t,y,MAP,Nstep,xNOL)ActualmodelwithPinfinity(t,y,MAP(:,k),Nstep,xNOL(k))
should only have two input arguments like this:
@(t,y)ActualmodelwithPinfinity(t,y,MAP(:,k),Nstep,xNOL(k))
You will then get other errors due to incompatible dimensions: it appears that you have used matrix operations everywhere in your code, but possibly you should be using array operations:
0 Comments
See Also
Categories
Find more on Probability Density Functions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!