-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsolvediscontODE_RT_Bev_Pem_eval_withPre.m
149 lines (119 loc) · 3.67 KB
/
solvediscontODE_RT_Bev_Pem_eval_withPre.m
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
function [ypred,yout,tout,v_end,t_vend,gammaout] = solvediscontODE_RT_Bev_Pem_eval_withPre(pars,Init,xdata,ydata)
global time_firstRT n_subpops useOnlyS
t_final_plot = 900*24; % Total time frame simulated
rt_on = [0,1,2,3,4]*24+time_firstRT; % in h - time of RT delivery
%odeset('Events',@eventsPat,'OutputFcn',@odeplot,'OutputSel',1,...
%'Refine',refine);
options = [];
% surviving fraction
if useOnlyS
S = pars(1:n_subpops);
else
error('Option useOnlyS=false is not supported!')
end
% Pretreatment
if any(xdata<0)
lambda = pars(4*n_subpops:5*n_subpops-1);
tout = [min(xdata)-14:1:0];
yout = zeros(3,length(tout));
yout(n_subpops+1:2*n_subpops,:) = Init(n_subpops+1:2*n_subpops)*exp(-(tout(1)-tout)*lambda);
Init(n_subpops+1:2*n_subpops)=yout(n_subpops+1:2*n_subpops,end);
else
tout = [0];
yout = [Init'];
end
% First fraction at time 0
v_viab = Init(n_subpops+1:2*n_subpops);
v_dying = Init(1:n_subpops);
v_dying = v_dying+(1-S).*v_viab;
v_viab = v_viab-(1-S).*v_viab;
gamma = Init(2*n_subpops+1);
Init = [v_dying,v_viab,gamma];
rt_on = rt_on(2:end);
tstart = 0;
tfinal = rt_on(1);
numfractions = length(rt_on);
for i = 1:numfractions+1
[t,y] = ode45(@model_ODE_pat_RT_Bev_Pem,[tstart:1:tfinal],Init,options,pars);
y = transpose(y);
t = transpose(t);
% Accumulate output and set the new initial conditions
v_viab = y(n_subpops+1:2*n_subpops,end)';
v_dying =y(1:n_subpops,end)';
gamma = y(1+2*n_subpops,end)';
v_dying = v_dying+(1-S).*v_viab;
v_viab = v_viab-(1-S).*v_viab;
Init=[v_dying,v_viab,gamma];
tout = [tout t(2:end)];
yout = [yout y(:,2:end)];
% Reassign start/end times
if i<numfractions
tstart = rt_on(i);
tfinal = rt_on(i+1);
else
tstart = tfinal;
tfinal = t_final_plot;
end
end
yall = sum(yout(1:2*n_subpops,:),1);
% Simulate more time if needed to show the full range until time to
% regrowth
countex = 0;
while max(yall)<1.1*yall(1) && countex <10
v_viab = y(n_subpops+1:2*n_subpops,end)';
v_dying =y(1:n_subpops,end)';
gamma = y(1+2*n_subpops,end)';
Init=[v_dying,v_viab,gamma];
tstart = tout(end);
tfinal = tstart+6*6*7*24;
sol = ode45(@model_ODE_pat_RT_Bev_Pem,[tstart:1:tfinal],Init,options,pars);
[t,y] = ode45(@model_ODE_pat_RT_Bev_Pem,[tstart:1:tfinal],Init,options,pars);
y = transpose(y);
t = transpose(t);
% Accumulate output and set the new initial conditions
v_viab = y(n_subpops+1:2*n_subpops,end)';
v_dying =y(1:n_subpops,end)';
gamma = y(1+2*n_subpops,end)';
Init=[v_dying,v_viab,gamma];
tout = [tout t(2:end)];
yout = [yout y(:,2:end)];
yall = nansum(yout(1:2*n_subpops,:),1);
countex = countex+1;
end
% Only use values which are not nan or inf
inds = ~isinf(yall)&~isnan(yall);
yall = yall(inds);
gammaout = yout(2*n_subpops+1,:);
yout = yout(1:2*n_subpops,inds);
tout = tout(inds);
[uni_tout,uniqueinds]=unique(tout);
if length(uniqueinds)~=length(yall)
yall=yall(uniqueinds);
tout = uni_tout;
yout = yout(1:2*n_subpops,uniqueinds);
gammaout =gammaout(inds);
gammaout =gammaout(uniqueinds);
end
% Interpolate at data points
try
ypred = interp1(tout,yall,xdata);
catch
disp(tout(inds));
end
% Get regrowth times - this has to be after the minimum vol has been
[~,minind] = min(yall);
v0_RT = yall(tout==0);
if ydata(end)>v0_RT
try
v_end = interp1(tout(minind:end),yall(minind:end),xdata(end));
catch
v_end = ydata(end);
end
else
v_end = 1.1*v0_RT;
end
try
t_vend = interp1(yall(minind:end),tout(minind:end),v_end);
catch
t_vend = 100000;
end