Skip to content
Snippets Groups Projects
Commit 4497cdf2 authored by Elin Nyman's avatar Elin Nyman
Browse files

Update code

parent 91429cd6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 142 additions and 1083 deletions
HypothesisA/LoadData.m
+ 4
− 11
View file @ 4497cdf2
......@@ -16,22 +16,21 @@ c4 = [13.15878 14.14623 18.20223]; l100_4 = [85.83585 89.46417 103.6841 105.3329
c6 = [14.30888 14.67227 18.2173]; l100_6 = [158.4887 154.0828 162.2422 255.6139 180.5526 190.212]; l1000_6 = [206.7315 213.3036 206.4169 348.8194 342.8387 306.1818];
c24 = [27.62139 34.96253 26.30299 57.61628 62.45974 69.22889]; l100_24 = [166.9165 147.1972 173.7321 630.5525 439.8181 446.5341]; l1000_24 = [212.4641 226.8367 232.186 558.759 568.0421 504.5358];
% Timepoints in the time response plot (figure 10)
% Timepoints in the time response plot
EXPDATA.timepoints = [0, 1, 2, 4, 6, 24];
% TNFA-response of the control in the time-response plot (figure 10)
% TNFA-response of the control in the time-response plot
EXPDATA.tnfaControl = [0, mean(c1), mean(c2), mean(c4), mean(c6), mean(c24)];
EXPDATA.semControl = [0, std(c1)/sqrt(3)/corrs3, std(c2)/sqrt(3)/corrs3, std(c4)/sqrt(3)/corrs3, std(c6)/sqrt(3)/corrs3, std(c24)/sqrt(6)/corrs6];
% TNFA-response of LPS conc. 100 ng/ml in the time-response plot (figure 10)
% TNFA-response of LPS conc. 100 ng/ml in the time-response plot
EXPDATA.tnfa100 = [0, mean(l100_1), mean(l100_2), mean(l100_4), mean(l100_6), mean(l100_24)];
EXPDATA.sem100 = [0, std(l100_1)/sqrt(3)/corrs3, std(l100_2)/sqrt(4)/corrs4, std(l100_4)/sqrt(6)/corrs6, std(l100_6)/sqrt(6)/corrs6, std(l100_24)/sqrt(6)/corrs6];
% TNFA-response of LPS conc. 1000 ng/ml in the time-response plot (figure 10)
% TNFA-response of LPS conc. 1000 ng/ml in the time-response plot
EXPDATA.tnfa1000 = [0, mean(l1000_1), mean(l1000_2), mean(l1000_4), mean(l1000_6), mean(l1000_24)];
EXPDATA.sem1000 = [0, std(l1000_1)/sqrt(3)/corrs3, std(l1000_2)/sqrt(6)/corrs6, std(l1000_4)/sqrt(6)/corrs6, std(l1000_6)/sqrt(6)/corrs6, std(l1000_24)/sqrt(6)/corrs6];
%x-axeln i figur 9 i Maria Linds master
EXPDATA.lps = [10,50,100,250,500,1000];
%LPS (ng/mL)
......@@ -42,18 +41,12 @@ c250 = [468.6055 382.417];
c500 = [465.6814 458.3753];
c1000 = [490.0401 485.2569];
%y-axeln i figur 9 i Maria Linds master
EXPDATA.tnfa = [mean(c10),mean(c50),mean(c100),mean(c250),mean(c500),mean(c1000)];
%residuals i figur 9 i Maria Linds master.
EXPDATA.sem = [std(c10)/sqrt(2)/corrs2,std(c50)/sqrt(2)/corrs2,std(c100)/sqrt(2)/corrs2,std(c250)/sqrt(2)/corrs2,std(c500)/sqrt(2)/corrs2,std(c1000)/sqrt(2)/corrs2];
%% Special SEM calculations
EXPDATA.sem100 = [0,mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),EXPDATA.sem100(end)]; %use SEM from dose-response (when higher)
%EXPDATA.sem100 = [0,EXPDATA.sem(3),EXPDATA.sem(3),EXPDATA.sem(3),EXPDATA.sem(3),EXPDATA.sem100(end)]; %use SEM from dose-response (when higher)
EXPDATA.sem1000 = [0,mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),EXPDATA.sem1000(end)]; %even out SEM
%EXPDATA.sem1000 = [0,EXPDATA.sem1000(end-1),EXPDATA.sem1000(end-1),EXPDATA.sem1000(end-1),EXPDATA.sem1000(end-1),EXPDATA.sem1000(end)]; %even out SEM
EXPDATA.semControl = [0,max(EXPDATA.semControl),max(EXPDATA.semControl),max(EXPDATA.semControl),max(EXPDATA.semControl),max(EXPDATA.semControl)]; %use SEM from dose-respose
%EXPDATA.semControl = [0,EXPDATA.sem(1),EXPDATA.sem(1),EXPDATA.sem(1),EXPDATA.sem(1),EXPDATA.sem(1)]; %use SEM from dose-respose
EXPDATA.sem = [mean(EXPDATA.sem), EXPDATA.sem(2), EXPDATA.sem(3), EXPDATA.sem(4), mean(EXPDATA.sem), mean(EXPDATA.sem)];
%EXPDATA.sem(end) = EXPDATA.sem1000(end); %use SEM from time-series
\ No newline at end of file
HypothesisA/LoadDrugData.m
+ 14
− 23
View file @ 4497cdf2
%% DEXAMETHASONE DATA (INHIBITION OF TNF-alfa DATA). MARIA LINDH %%
%% DEXAMETHASONE DATA, MARIA LINDH %%
DRUGDATA = [];
% Timepoints in the drug response time-series (figure 12)
% Timepoints in the drug response time-series
DRUGDATA.timepoints = [0, 1, 2, 4, 6, 24];
%%% https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation
......@@ -22,25 +22,23 @@ c4 = [62.75599 57.82652]; l4 = [14.40909 17.0369 8.907761]; h4 = [2.258175 5.61
c6 = [194.3419 172.892]; l6 = [39.0532 43.15816 27.96148]; h6 = [16.71108 21.23619 12.19983]; l6p = [122.1991 197.0661 117.1118]; h6p = [50.12288 34.81426 32.16351];
c24 = [450.8582 485.3612]; l24 = [55.49924 53.36025 45.31882]; h24 = [29.17651 26.14517 24.43495]; l24p = [257.5461 391.8945 307.7493]; h24p = [96.68551 84.14325 97.12671];
% TNFA-response of the positive control (no drug) in the time-series (figure 12)
% TNFA-response of the positive control (no drug) in the time-series
DRUGDATA.tnfa100Control = [0, 0, mean(c2), mean(c4), mean(c6), mean(c24)];
DRUGDATA.sem100Control = [0, 0, std(c2)/sqrt(2)/corrs2, std(c4)/sqrt(2)/corrs2, std(c6)/sqrt(2)/corrs2, std(c24)/sqrt(2)/corrs2];
% TNFA-response of 0.3 �M dexamethasone in the time-series plot (figure 12)
% TNFA-response of 0.3 uM dexamethasone in the time-series plot
DRUGDATA.tnfaLowdex = [0, 0, 0, mean(l4), mean(l6), mean(l24)];
DRUGDATA.semLowdex = [0, 0, 0, std(l4)/sqrt(3)/corrs3, std(l6)/sqrt(3)/corrs3, std(l24)/sqrt(3)/corrs3];
% TNFA-response of 3 �M dexamethasone in the time-series plot (figure 12)
% TNFA-response of 3 uM dexamethasone in the time-series plot
DRUGDATA.tnfaHighdex = [0, 0, 0, mean(h4), mean(h6), mean(h24)];
DRUGDATA.semHighdex = [0, 0, 0, std(h4)/sqrt(2)/corrs2, std(h6)/sqrt(3)/corrs3, std(h24)/sqrt(3)/corrs3];
% TNFA-response of 0.3 �M dexamethasone (pre-treatment) in the time-series.
% (fig12)
% TNFA-response of 0.3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.tnfaLowdex_pre = [0, 0, 0, mean(l4p), mean(l6p), mean(l24p)];
DRUGDATA.semLowdex_pre = [0, 0, 0, std(l4p)/sqrt(3)/corrs3, std(l6p)/sqrt(3)/corrs3, std(l24p)/sqrt(3)/corrs3];
% TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig12)
% TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.tnfaHighdex_pre = [0, 0, 0, mean(h4p), mean(h6p), mean(h24p)];
DRUGDATA.semHighdex_pre = [0, 0, 0, std(h4p)/sqrt(3)/corrs3, std(h6p)/sqrt(3)/corrs3, std(h24p)/sqrt(3)/corrs3];
......@@ -51,7 +49,6 @@ DRUGDATA.semHighdex = [0,0,0,max(DRUGDATA.semHighdex),max(DRUGDATA.semHighdex),m
DRUGDATA.semLowdex_pre = [0,0,0,max(DRUGDATA.semLowdex_pre),max(DRUGDATA.semLowdex_pre),max(DRUGDATA.semLowdex_pre)];
DRUGDATA.semHighdex_pre = [0,0,0,max(DRUGDATA.semHighdex_pre),max(DRUGDATA.semHighdex_pre),max(DRUGDATA.semHighdex_pre)];
%% Fig 13 data (pre-treatment at different timepoints)
%Time of LPS incubation (hours) control control control pre-treatment 1 pre-treatment 1 pre-treatment 1 pre-treatment 2 pre-treatment 2 pre-treatment 2 pre-treatment 3 pre-treatment 3 pre-treatment 3
c0 = [2.47024 1.00613]; p1_0 = [0]; p2_0 = [0]; p3_0 = [0.656634];
c1 = [9.552152 9.736456 5.547666]; p1_1 = [2.686642]; p2_1 = [0]; p3_1 = [0];
......@@ -60,36 +57,30 @@ c4 = [103.4877 172.0678 93.0803]; p1_4 = [65.6666 54.6981 55.74894]; p2_4 = [54.
c6 = [152.6309 220.9061 122.7121]; p1_6 = [89.44148 67.00409 85.37128]; p2_6 = [77.26416 86.26146 125.5514]; p3_6 = [32.67914 37.54742 61.43369];
%Annova
Matr = [c6; p3_6; p2_6; p1_6]';
[~,~,stats] = anova1(Matr);
[c,~,~,gnames] = multcompare(stats,'CType','bonferroni');
c
%Matr = [c6; p3_6; p2_6; p1_6]';
%[~,~,stats] = anova1(Matr);
%[c,~,~,gnames] = multcompare(stats,'CType','bonferroni');
% TIMEPOINTS
DRUGDATA.timepoints2 = [0, 1, 2, 4, 6];
% CONTROL: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% CONTROL: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre_control_tnfa = [mean(c0), mean(c1), mean(c2), mean(c4), mean(c6)];
DRUGDATA.pre_control_sem = [std(c0)/sqrt(2)/corrs2, std(c1)/sqrt(3)/corrs3, std(c2)/sqrt(3)/corrs3, std(c4)/sqrt(3)/corrs3, std(c6)/sqrt(3)/corrs3];
% PRETREATMENT1: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT1: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre1_tnfa = [0, 0, mean(p1_2), mean(p1_4), mean(p1_6)];
DRUGDATA.pre1_sem = [0, 0, std(p1_2)/sqrt(3)/corrs3, std(p1_4)/sqrt(3)/corrs3, std(p1_6)/sqrt(3)/corrs3];
% PRETREATMENT2: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT2: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre2_tnfa = [0, 0, mean(p2_2), mean(p2_4), mean(p2_6)];
DRUGDATA.pre2_sem = [0, 0, std(p2_2)/sqrt(3)/corrs3, std(p2_4)/sqrt(3)/corrs3, std(p2_6)/sqrt(3)/corrs3];
% PRETREATMENT3: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT3: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre3_tnfa = [0, 0, mean(p3_2), mean(p3_4), mean(p3_6)];
DRUGDATA.pre3_sem = [0, 0, std(p3_2)/sqrt(3)/corrs3, std(p3_4)/sqrt(3)/corrs3, std(p3_6)/sqrt(3)/corrs3];
%Set SEM to MAX(SEM)
%DRUGDATA.pre_control_sem = [max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem)];
DRUGDATA.pre1_sem = [0,0,max(DRUGDATA.pre1_sem),max(DRUGDATA.pre1_sem),max(DRUGDATA.pre1_sem)];
DRUGDATA.pre2_sem = [0,0,max(DRUGDATA.pre2_sem),max(DRUGDATA.pre2_sem),max(DRUGDATA.pre2_sem)];
DRUGDATA.pre3_sem = [0,0,max(DRUGDATA.pre3_sem),max(DRUGDATA.pre3_sem),max(DRUGDATA.pre3_sem)];
HypothesisA/Macrophage_simple.txt
+ 1
− 1
View file @ 4497cdf2
********** MODEL NAME
Macrophage_simple
Hypothesis A
********** MODEL NOTES
time in hours, concentrations in uM
......
HypothesisA/Optimization_simple.m
+ 15
− 45
View file @ 4497cdf2
clear all
close all
global EXPDATA
......@@ -10,7 +9,6 @@ global TNF_Timepoints
global DRUGDATA
global TNF_Timepoints2
global FID
global prev_cost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD THE MODEL
......@@ -27,58 +25,30 @@ TNF_Timepoints = DRUGDATA.timepoints;
TNF_Timepoints2 = DRUGDATA.timepoints2;
FID = fopen('allGoodValues.dat','wt');
prev_cost = 300;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THE OPTIMIZATION
psOpts = optimoptions(@particleswarm,'Display','iter');
saOpts = optimoptions(@simulannealbnd, 'HybridFcn',@fmincon,'Display','iter');%'TolFun', 1e-6, 'StallIterLimit', 1000,
saOpts = optimoptions(@simulannealbnd, 'HybridFcn',@fmincon,'Display','iter');
X = log(paramValues(1:16))';
X = [-3.0891 2.2374 -1.6336 0.3039 -0.8985 -1.3299 4.4345 6.2120 2.4568 2.7877 3.5930 -2.5122 -1.0078 4.9256 -4.2607 2];
X = [-2.6181 2.1886 -1.9698 0.3015 1.3786 -0.5108 4.6628 6.7670 1.6537 3.1504 2.3178 -2.0874 -0.2298 4.4830 -6.0119 1.7323];
X = [-0.73342 4.7795 0.07908 -0.48928 2.7646 2.749 6.7876 5.9975 4.4351 0.95724 1.4047 -0.63016 -0.95093 4.976 -6.289 4.3447];
nParams = length(X);
%%
for i=1:20
i
%Define upper and lower bounds
lb = ones(16,1) * 1e-3;
%lb(3) = 0.19 * 0.5;
%lb(4) = 2.7 * 0.5;
%lb(5) = 0.5*60; %kon 30
%lb(6) = 0.0001*60; %koff 0.06
%lb(6) = 0.0034*60*60; 12
lb=log(lb);
ub = ones(16,1) * 1e3;
%ub(3) = 0.19 * 2;
%ub(4) = 2.7 * 2;
%ub(5) = 1*60; %kon 60
%ub(6) = 0.01 * 60; %koff 0.6
%ub(6) = 0.1; %!!!
%ub(6) = 0.007*60*60; %25
%ub(6) = 5
ub=log(ub);
format long
format compact
%Define upper and lower parameter boundaries
lb = ones(16,1) * 1e-3;
lb=log(lb);
ub = ones(16,1) * 1e3;
ub=log(ub);
format long
format compact
%[optParamPS, minfunPS]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
%[X, FVAL]=simulannealbnd(@cost_simple, optParamPS, lb,ub,saOpts);
[X, FVAL]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
%[X, FVAL]=simulannealbnd(@cost_simple, X, lb,ub,saOpts);
save(sprintf('opt(%.2f), %s.mat',FVAL, datestr(now,'yymmdd-HHMMSS')),'X')
end
i
[optParamPS, minfunPS]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
[X, FVAL]=simulannealbnd(@cost_simple, optParamPS, lb,ub,saOpts);
save(sprintf('opt(%.2f), %s.mat',FVAL, datestr(now,'yymmdd-HHMMSS')),'X')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fclose(FID)
HypothesisA/cost_simple.m
+ 30
− 55
View file @ 4497cdf2
......@@ -6,29 +6,23 @@ function [cost] = cost_simple(param, shouldIPlot)
global Model
global TNF_Timepoints
global TNF_Timepoints2
global model1_EXPDATA
global DRUGDATA
global FID
global prev_cost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SIMULATION
param = exp(param);
%param(end) = 0;
%%STEADY STATE SIMULATION (LPS = Dexa = 0)
param = [param, 0, 0];
param = exp(param);
%%STEADY STATE SIMULATION (LPS = Dexa = 0)
param = [param, 0, 0];
try
try
ss_simulation = IQMPsimulate(Model, 700, initCon, paramNames, param);
ss_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
ss_simulation = IQMPsimulate(Model, 700, initCon, paramNames, param);
ss_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%% LPS simulations from Maria Linds master thesis (figure 10 in her report).
%% LPS simulations
param(end) = 0;
param(end-1) = 0;
......@@ -48,37 +42,36 @@ LPS_1000_simulation = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.stateval
param(end-1) = 0;
doseresponseLPS = [LPS_10_simulation.variablevalues(end,2),LPS_50_simulation.variablevalues(end,2),LPS_100_simulation.variablevalues(end,2),LPS_250_simulation.variablevalues(end,2),LPS_500_simulation.variablevalues(end,2),LPS_1000_simulation.variablevalues(end,2)];
%% Dexamethasone simulations using data from Maria Linds master thesis (figure 12 in her report).
%% Dexamethasone simulations
% PRE-TREATMENT SIMULATIONS - macrophages are pre-treated with 0.3 or 3 �g
% PRE-TREATMENT SIMULATIONS - macrophages are pre-treated with 0.3 or 3 ug
% Dexa for 1 hour
param(end) = 0.3;
Dexa_simulation_pre300 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
param(end) = 3;
Dexa_simulation_pre3000 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
% 0 �M DEXAMETHASONE SIMULATION (positive control), only 100 ng/ml LPS
% 0 uM DEXAMETHASONE SIMULATION (positive control), only 100 ng/ml LPS
param = param(1:(end-2));
param = [param, 100, 0];
Dexa_simulation1 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end,:), paramNames, param);
% 0.3 �M DEXAMETHASONE SIMULATION
% 0.3 uM DEXAMETHASONE SIMULATION
param(end) = 0.3;
Dexa_simulation2 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end, :), paramNames, param);
% 3 �M DEXAMETHASONE SIMULATION
% 3 uM DEXAMETHASONE SIMULATION
param(end) = 3;
Dexa_simulation3 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end, :), paramNames, param);
% 0.3 �M DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0; % 1 % of 0.3 uM remains
% 0.3 uM DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0;
Dexa_simulation4 = IQMPsimulate(Model, TNF_Timepoints, Dexa_simulation_pre300.statevalues(end, :), paramNames, param);
% 3 �M DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0; % 1 % of 3 uM remains
% 3 uM DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0;
Dexa_simulation5 = IQMPsimulate(Model, TNF_Timepoints, Dexa_simulation_pre3000.statevalues(end, :), paramNames, param);
%% Dexamethasone simulations using data from Maria Linds master thesis (figure 13 in her report).
%Experiments with pre-treatment at different timepoints.
......@@ -87,14 +80,14 @@ param = param(1:(end-2));
param = [param, 100, 0];
Pre_control_simulation = IQMPsimulate(Model, TNF_Timepoints2, ss_simulation.statevalues(end,:), paramNames, param);
% % 1 hour stimulation with 3 �g dexa (pre-treatment)
% % 1 hour stimulation with 3 ug dexa (pre-treatment)
param = param(1:(end-2));
param = [param, 0, 3];
one_hour_simulation = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
one_hour_simulation.statevalues(end, end) = 0; %TNFa 0, buffer change
% % Pre-treatment 1
param(end) = 0; % 1 % of 3 uM remains
param(end) = 0;
twentytwo_hour_simulation = IQMPsimulate(Model, 22, one_hour_simulation.statevalues(end,:), paramNames, param);
param(end-1) = 100;
twentytwo_hour_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
......@@ -110,17 +103,10 @@ Pre2_simulation = IQMPsimulate(Model, TNF_Timepoints2, sixteen_hour_simulation.s
% % Pre-treatment 3
Pre3_simulation = IQMPsimulate(Model, TNF_Timepoints2, one_hour_simulation.statevalues(end,:), paramNames, param);
% % Pre-treatment long
%param(end-1) = 0;
%many_hour_simulation = IQMPsimulate(Model, 168, one_hour_simulation.statevalues(end,:), paramNames, param);
%param(end-1) = 100;
%many_hour_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%PreLong_simulation = IQMPsimulate(Model, TNF_Timepoints2, many_hour_simulation.statevalues(end,:), paramNames, param);
catch myError
cost = 1e6;
return
end
catch myError
cost = 1e6;
return
end
Control_simval = Pre_control_simulation.variablevalues(:,2);
Pre1_simval = Pre1_simulation.variablevalues(:,2);
......@@ -134,21 +120,15 @@ Drug_simvalues4 = Dexa_simulation4.variablevalues(:,2);
Drug_simvalues5 = Dexa_simulation5.variablevalues(:,2);
%scale day-to-day differences between sets of data linearly
%scale=lscov([Control_simval; Pre1_simval; Pre2_simval; Pre3_simval], [DRUGDATA.pre_control_tnfa'; DRUGDATA.pre1_tnfa'; DRUGDATA.pre2_tnfa'; DRUGDATA.pre3_tnfa']);
scale=lscov([Drug_simvalues1; Drug_simvalues2; Drug_simvalues3; Drug_simvalues4; Drug_simvalues5], [DRUGDATA.tnfa100Control'; DRUGDATA.tnfaLowdex'; DRUGDATA.tnfaHighdex'; DRUGDATA.tnfaLowdex_pre'; DRUGDATA.tnfaHighdex_pre']);
%scale_exp3=lscov([LPS_0_simulation.variablevalues(:,2); LPS_100_simulation.variablevalues(:,2); LPS_1000_simulation.variablevalues(:,2)], [EXPDATA.tnfaControl'; EXPDATA.tnfa100'; EXPDATA.tnfa1000']);
cost1 = sum((DRUGDATA.pre_control_tnfa - scale.*Control_simval').^2 ./ (DRUGDATA.pre_control_sem).^2);
cost2 = sum((DRUGDATA.pre1_tnfa(3:end) - scale.*Pre1_simval(3:end)').^2 ./ (DRUGDATA.pre1_sem(3:end)).^2);
cost3 = sum((DRUGDATA.pre2_tnfa(3:end) - scale.*Pre2_simval(3:end)').^2 ./ (DRUGDATA.pre2_sem(3:end)).^2);
cost4 = sum((DRUGDATA.pre3_tnfa(3:end) - scale.*Pre3_simval(3:end)').^2 ./ (DRUGDATA.pre3_sem(3:end)).^2);
%if abs(DRUGDATA.pre_control_tnfa(end) - scale.*PreLong_simulation.variablevalues(end,2)) > 20
% costextra = 50;
%else
% costextra = 0;
%end
% Not used in training
%cost1 = sum((DRUGDATA.pre_control_tnfa - scale.*Control_simval').^2 ./ (DRUGDATA.pre_control_sem).^2);
%cost2 = sum((DRUGDATA.pre1_tnfa(3:end) - scale.*Pre1_simval(3:end)').^2 ./ (DRUGDATA.pre1_sem(3:end)).^2);
%cost3 = sum((DRUGDATA.pre2_tnfa(3:end) - scale.*Pre2_simval(3:end)').^2 ./ (DRUGDATA.pre2_sem(3:end)).^2);
%cost4 = sum((DRUGDATA.pre3_tnfa(3:end) - scale.*Pre3_simval(3:end)').^2 ./ (DRUGDATA.pre3_sem(3:end)).^2);
% Used in training
cost5 = sum((DRUGDATA.tnfa100Control(3:end) - scale.*Drug_simvalues1(3:end)').^2 ./ (DRUGDATA.sem100Control(3:end)).^2);
cost6 = sum((DRUGDATA.tnfaLowdex(4:end) - scale.*Drug_simvalues2(4:end)').^2 ./ (DRUGDATA.semLowdex(4:end)).^2);
cost7 = sum((DRUGDATA.tnfaHighdex(4:end) - scale.*Drug_simvalues3(4:end)').^2 ./ (DRUGDATA.semHighdex(4:end)).^2);
......@@ -161,19 +141,14 @@ costLPS1000 = sum((EXPDATA.tnfa1000(2:end) - LPS_1000_simulation.variablevalues(
costLPSdose = sum((EXPDATA.tnfa - doseresponseLPS).^2 ./ (EXPDATA.sem).^2);
%size([DRUGDATA.pre_control_tnfa,DRUGDATA.pre1_tnfa(3:end),DRUGDATA.pre2_tnfa(3:end),DRUGDATA.pre3_tnfa(3:end),DRUGDATA.tnfa100Control(3:end),DRUGDATA.tnfaLowdex(4:end),DRUGDATA.tnfaHighdex(4:end),DRUGDATA.tnfaHighdex_pre(4:end),EXPDATA.tnfa,EXPDATA.tnfa1000(2:end),EXPDATA.tnfa100(2:end),EXPDATA.tnfaControl(2:end),DRUGDATA.tnfaLowdex_pre(4:end)],2)
%costTNF = cost1 + cost2 + cost3 + cost4 + cost5 + cost6 + cost7 + cost8 + cost9;
costTNF = cost5 + cost6 + cost7 + cost8 + cost9;
costLPS = costLPS0 + costLPS100 + costLPS1000;
cost = costTNF + costLPSdose + costLPS; %+ costextra
%cost = cost1 + cost2 + cost3 + cost4;
%cost1,cost2,cost3,cost4,cost5,cost6,cost7,cost8,cost9,costLPS0,costLPS100,costLPS1000,costLPSdose
%% CHI-SQUARE TEST
if cost < 65
%if cost < prev_cost
prev_cost = cost;
fprintf(FID,'%4.10f %10.10f ',[cost, param(1:end-2), scale]); fprintf(FID,'\n');
end
......
HypothesisA/plot_many_params.m
+ 11
− 165
View file @ 4497cdf2
%plot_many_params
close all
%% LOAD MODEL AND DATA
Model = 'Macrophage_simple';
optModel = IQMmodel(strcat(Model,'.txt'));
......@@ -33,17 +33,13 @@ dexa03 = [213,94,0]./256; %red
dexa3_pre3 = [86,180,233]./256; %sky blue
dexa3 = [0,114,178]./256; %blue
dexa3_pre2 = [0,158,115]./256; %green
%dexa3_pre1 = [60,218,175]./256; %light green
%dexa3_pre1 = [80,238,195]./256; %light green
dexa3_pre1 = [80,0,0]./256; %brown
thres = 69; %all data chi2inv(0.95,51)
thres = 65; %validation pretreatments chi2inv(0.95,48)
thres = 52; %validation pretreatments chi2inv(0.95,37)
thres = 52;
%% LOAD PARAMS
values1=load('allGoodValues1.dat'); %no exp D, constrained koff
values2=load('allGoodValues2.dat'); %no exp D, no constrain
values3=load('allGoodValues3.dat'); %no exp D, no constrain
values4=load('allGoodValues4.dat'); %no exp D, no constrain
......@@ -51,12 +47,11 @@ values5=load('allGoodValues5.dat'); %no exp D, no constrain
values6=load('allGoodValues6.dat'); %no exp D, no constrain
values = unique([values6; values5; values4; values3; values2],'rows');
%values = unique([values2],'rows');
ind=find(values(:,1)<thres);
res = 100;
res = 1000;
all_param = ones(round(size(ind,1)/res), size(values1,2));
all_param = ones(round(size(ind,1)/res), size(values2,2));
size(all_param)
l=1;
......@@ -138,37 +133,6 @@ for i = 1:res:size(ind, 1)
X = X(1:(end-2));
X = [X, 100, 0];
pre3_sim = IQMPsimulate(Model, 7, one_hour_sim.statevalues(end,:), paramNames, X);
if 1.07*pre1_sim.variablevalues(end,2) < Dexa_simulation.variablevalues(round(1001/25*7),2)
best_pre1 = scale.*pre1_sim.variablevalues(:,2);
best_pre2 = scale.*pre2_sim.variablevalues(:,2);
best_pre3 = scale.*pre3_sim.variablevalues(:,2);
best_control = scale.*Dexa_simulation.variablevalues(:,2);
best_10 = LPS_10_simulation.variablevalues(end,2);
best_50 = LPS_50_simulation.variablevalues(end,2);
best_100 = LPS_100_simulation.variablevalues(end,2);
best_250 = LPS_250_simulation.variablevalues(end,2);
best_500 = LPS_500_simulation.variablevalues(end,2);
best_1000 = LPS_1000_simulation.variablevalues(end,2);
best_0 = LPS_0_simulation.variablevalues(:,2);
best_100_x = LPS_100_simulation.variablevalues(:,2);
best_1000_x = LPS_1000_simulation.variablevalues(:,2);
best_Dexa = scale.*Dexa_simulation.variablevalues(:,2);
best_sim2 = scale.*Model_simulation2.variablevalues(:,2);
best_sim3 = scale.*Model_simulation3.variablevalues(:,2);
best_sim4 = scale.*Model_simulation4.variablevalues(:,2);
best_sim5 = scale.*Model_simulation5.variablevalues(:,2);
%IKB
basal = 0;%one_hour_sim.statevalues(1,5);
maxi = one_hour_sim.statevalues(end,5);
best_IKB0 = (one_hour_sim.statevalues(:,5)-basal)./maxi;
best_IKBbrown = (LPS_100_simulation.statevalues(:,5)-basal)./maxi;
best_IKBpink = (pre3_sim.statevalues(:,5)-basal)./maxi;
best_IKBgreen = (Model_simulation3.statevalues(:,5)-basal)./maxi;
end
if l == 1
......@@ -227,7 +191,7 @@ for i = 1:res:size(ind, 1)
p3n_max = (pre3_sim.statevalues(:,8)-basal)./maxi;
%IKB
basal = 0;%one_hour_sim.statevalues(1,5);
basal = 0;
maxi = one_hour_sim.statevalues(end,5);
IKB0_max = (one_hour_sim.statevalues(:,5)-basal)./maxi;
IKB0_min = (one_hour_sim.statevalues(:,5)-basal)./maxi;
......@@ -311,7 +275,7 @@ for i = 1:res:size(ind, 1)
end
%IKB
basal = 0;%one_hour_sim.statevalues(1,5);
basal = 0;
maxi = one_hour_sim.statevalues(end,5);
for j = 1:size(IKB0_min, 1)
......@@ -338,43 +302,9 @@ for i = 1:res:size(ind, 1)
IKBgreen_max(j) = max(IKBgreen_max(j), (Model_simulation3.statevalues(j,5)-basal)./maxi);
IKBgreen_min(j) = min(IKBgreen_min(j), (Model_simulation3.statevalues(j,5)-basal)./maxi);
end
%figure(105)
%hold on
%plot(0:0.001:1, (one_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%plot(22:7/1000:29, (pre1_sim.statevalues(:,9)-basal)./maxi, 'Color', dexa3_pre1, 'LineWidth', 0.5);
%plot(16:7/1000:23, (pre2_sim.statevalues(:,9)-basal)./maxi, 'Color',dexa3_pre2, 'LineWidth', 0.5);
%plot(1:7/1000:8, (pre3_sim.statevalues(:,9)-basal)./maxi, 'Color', dexa3_pre3, 'LineWidth', 0.5);
%plot(1:15/1000:16, (sixteen_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%plot(1:21/1000:22, (twentytwo_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%figure(101)
%for i = 1:7
% basal = one_hour_sim.statevalues(1,i+5);
% %subplot(4,2,i), plot(LPS_100_simulation.time, LPS_100_simulation.statevalues(:,i+5), 'Color', LPS100, 'LineWidth', 0.5);
% hold on
% subplot(4,2,i),plot(1:7/1000:8, pre3_sim.statevalues(:,i+5)-basal, 'Color', dexa3_pre3, 'LineWidth', 0.5);
% subplot(4,2,i),plot(16:7/1000:23, pre2_sim.statevalues(:,i+5)-basal, 'Color',dexa3_pre2, 'LineWidth', 0.5);
% subplot(4,2,i),plot(22:7/1000:29, pre1_sim.statevalues(:,i+5)-basal, 'Color', dexa3_pre1, 'LineWidth', 0.5);
% %subplot(4,2,i),plot(1:25/1000:26, Model_simulation3.statevalues(:,i+5), 'Color', dexa3, 'LineWidth', 0.5);
% subplot(4,2,i), plot(0:0.001:1, one_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
% subplot(4,2,i),plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
% subplot(4,2,i), plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
%end
end
% figure(3)
% hold on
% plot(Dexa_simulation.time, scale_exp1.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 0.5);
% plot(pre1_sim.time, scale_exp1.*pre1_sim.variablevalues(:,2), 'Color', dexa3_pre1, 'LineWidth', 0.5);
% plot(pre2_sim.time, scale_exp1.*pre2_sim.variablevalues(:,2), 'Color', dexa3_pre2, 'LineWidth', 0.5);
% plot(pre3_sim.time, scale_exp1.*pre3_sim.variablevalues(:,2), 'Color', dexa3_pre3, 'LineWidth', 0.5);
end
%end
%% PLOT
figure(1)
subplot(2,2,1)
......@@ -383,8 +313,6 @@ errorbar(EXPDATA.lps, EXPDATA.tnfa, EXPDATA.sem, 'ks','LineWidth', 2.5, 'MarkerF
h = fill([EXPDATA.lps,flip(EXPDATA.lps)], [LPS_10_min,LPS_50_min,LPS_100_min,LPS_250_min,LPS_500_min,LPS_1000_min,LPS_1000_max,LPS_500_max,LPS_250_max,LPS_100_max,LPS_50_max,LPS_10_max], 'k');
set(h,'facealpha',.6,'EdgeColor','none')
%plot(EXPDATA.lps, [best_10,best_50,best_100,best_250,best_500,best_1000], 'Color', LPS0, 'LineWidth', 2);
title('A) LPS dose-response');
xlabel('LPS (ng/ml)');
ylabel('TNF (pg/mg tissue)');
......@@ -404,10 +332,6 @@ set(h,'facealpha',.6,'EdgeColor','none')
h = fill([LPS_0_simulation.time,flip(LPS_0_simulation.time)], ([LPS_0_min;flip(LPS_0_max)]'), LPS0);
set(h,'facealpha',.6,'EdgeColor','none')
%plot(LPS_0_simulation.time, best_0, 'Color', LPS0, 'LineWidth', 2);
%plot(LPS_100_simulation.time, best_100_x, 'Color', LPS100, 'LineWidth', 2);
%plot(LPS_1000_simulation.time, best_1000_x, 'Color', LPS1000, 'LineWidth', 2);
axis([0, 25, 0, 640]);
title('B) LPS time-response');
xlabel('Time (hours)');
......@@ -434,12 +358,6 @@ set(h,'facealpha',.6,'EdgeColor','none')
h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([sim5_min;flip(sim5_max)]'), dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
%plot(Dexa_simulation.time, best_Dexa, 'Color', LPS100, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim2, 'Color', dexa03, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim3, 'Color', dexa3, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim4, 'Color', dexa03_pre, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim5, 'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 25, 0, 640]);
title('C) Dexamethasone + LPS');
xlabel('Time (hours)');
......@@ -448,7 +366,6 @@ set(gca,'FontSize',14);
LEG = legend('control (only LPS)', '0.3 uM (removed)', '3 uM (removed)', '0.3 uM', '3 uM', 'Location', 'northwest');
set(LEG,'FontSize',14);
subplot(2,2,4)
hold on
h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([IKBgreen_min;flip(IKBgreen_max)]'), dexa3);
......@@ -459,15 +376,8 @@ h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([IKBbrown_min;flip(
set(h,'facealpha',.5,'EdgeColor','none')
h = fill([-1:0.001:0,flip(-1:0.001:0)], ([IKB0_min;flip(IKB0_max)]'), dexa3);
set(h,'facealpha',.5,'EdgeColor','none')
%plot(-1:0.001:0, best_IKB0, 'Color', dexa3, 'LineWidth', 2);
%plot(-1:0.001:0, best_IKB0, 'Color', dexa3, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_IKBbrown, 'Color', LPS100, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_IKBgreen, 'Color', dexa3, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_IKBpink, 'Color', dexa3_pre3, 'LineWidth', 2);
plot([0 0], [0 20], 'Color', LPS0, 'LineWidth', 0.5);
%xlim([-1 5])
axis([-1, 5, 0, IKBgreen_max(201)]);
title('D) Mechanism of sustained response');
xlabel('Time (hours)');
......@@ -477,81 +387,18 @@ LEG = legend('3 uM Dexa', '3 uM Dexa (removed)', 'control (LPS)', 'Location', 'n
set(LEG,'FontSize',14);
figure(2)
subplot(1,2,1)
hold on
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre1_min;flip(pre1_max)]'), dexa3_pre1);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre2_min;flip(pre2_max)]'), dexa3_pre2);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([Dexa_min;flip(Dexa_max)]'), LPS100);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre3_min;flip(pre3_max)]'), dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
plot(pre1_sim.time, best_pre1, 'Color', dexa3_pre1, 'LineWidth', 2);
plot(pre1_sim.time, best_pre2, 'Color', dexa3_pre2, 'LineWidth', 2);
plot(pre1_sim.time, best_pre3, 'Color', dexa3_pre3, 'LineWidth', 2);
plot(Dexa_simulation.time, best_control, 'Color', LPS100, 'LineWidth', 2);
axis([0, 7, 0, 260]);
title('Model predictions: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
set(gca,'FontSize',14);
%LEG = legend('control (only LPS)', '22 h wash before LPS', '16 h wash before LPS', 'dexa removed before LPS', 'Location', 'northwest');
%set(LEG,'FontSize',14);
subplot(1,2,2)
hold on
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 2.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 8);
axis([0, 7, 0, 260]);
title('Experimental data: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
set(gca,'FontSize',14);
LEG = legend('control (only LPS)', '22 h wash before LPS', '16 h wash before LPS', 'dexa removed before LPS', 'Location', 'northwest');
set(LEG,'FontSize',14);
%figure(101)
%for i = 1:7
% subplot(4,2,i),title(statenames(i+5))
%end
%scale2=lscov(mean([Dexa_min(DRUGDATA.timepoints2*40 + 1),Dexa_max(DRUGDATA.timepoints2*40 + 1)],2),DRUGDATA.pre_control_tnfa')
%scale2=0.84447
scale2=1
figure(102)
subplot(2,1,1)
hold on
h = fill([0:0.025:6, flip(0:0.025:6)], scale2.*[Dexa_min(1:241);flip(Dexa_max(1:241))]', LPS100);
h = fill([0:0.025:6, flip(0:0.025:6)], [Dexa_min(1:241);flip(Dexa_max(1:241))]', LPS100);
set(h,'facealpha',.6,'EdgeColor','none')
%h=fill([0:0.001:1, flip(0:0.001:1)],[one_min; flip(one_max)]', LPS0);
%set(h,'facealpha',.6,'EdgeColor','none')
h=fill([1:7/1000:8, flip(1:7/1000:8)],scale2.*[pre3_min; flip(pre3_max)]', dexa3_pre3);
h=fill([1:7/1000:8, flip(1:7/1000:8)],[pre3_min; flip(pre3_max)]', dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
h=fill([17:7/1000:24, flip(17:7/1000:24)],scale2.*[pre2_min; flip(pre2_max)]', dexa3_pre2);
h=fill([17:7/1000:24, flip(17:7/1000:24)],[pre2_min; flip(pre2_max)]', dexa3_pre2);
set(h,'facealpha',.6,'EdgeColor','none')
%h=fill([1:22/1000:23, flip(1:22/1000:23)],[twentytwo_min; flip(twentytwo_max)]', LPS0);
%set(h,'facealpha',.6,'EdgeColor','none')
h=fill([23:7/1000:30, flip(23:7/1000:30)],scale2.*[pre1_min; flip(pre1_max)]', dexa3_pre1);
h=fill([23:7/1000:30, flip(23:7/1000:30)],[pre1_min; flip(pre1_max)]', dexa3_pre1);
set(h,'facealpha',.6,'EdgeColor','none')
%errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 2.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 8);
%errorbar(DRUGDATA.timepoints2 + 1, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 8);
%errorbar(DRUGDATA.timepoints2 + 17, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 8);
%errorbar(DRUGDATA.timepoints2 + 23, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 8);
axis([0, 30, 0, max([Dexa_max(241),pre1_max(end)])+1]);
% hold on
% plot(0:0.001:1, (one_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
% plot(22:7/1000:29, (pre1_sim.statevalues(:,9)-basal)/maxi, 'Color', dexa3_pre1, 'LineWidth', 0.5);
% plot(16:7/1000:23, (pre2_sim.statevalues(:,9)-basal)/maxi, 'Color',dexa3_pre2, 'LineWidth', 0.5);
% plot(1:7/1000:8, (pre3_sim.statevalues(:,9)-basal)/maxi, 'Color', dexa3_pre3, 'LineWidth', 0.5);
%plot(1:15/1000:16, (sixteen_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
% plot(1:21/1000:22, (twentytwo_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
title('Model predictions: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
......@@ -574,6 +421,7 @@ plot([29 29], [210-5 210+5], 'Color', LPS0, 'LineWidth', 0.5);
plot([6 6], [120-5 120+5], 'Color', LPS0, 'LineWidth', 0.5);
text(6.5, 100, 'p = 0.0047','FontSize',14)
text(27, 190, 'p = 0.034','FontSize',14)
title('Experimental data: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
......@@ -582,10 +430,8 @@ LEG = legend('control (only LPS)', 'dexa removed before LPS', '16 h wash before
set(LEG,'FontSize',14);
orderp = [1,3,4,8,7,2,11,10,9,14,15,12,13,5,6,16,17];
figure(103)
figure(3)
semilogy(1:size(all_param,2)-1,all_param(:,1+orderp),'-o')
%title('Acceptable parameter values');
%xlabel('Name of parameter');
set(gca,'XTick',1:size(all_param,2)-1)
set(gca,'XTickLabel',[paramNames(orderp(1:end-1));'scale'])
set(gca,'XTickLabelRotation',45)
......
HypothesisA/simple_plot.m deleted 100644 → 0
+ 0
− 226
View file @ 91429cd6
function [ a ] = simple_plot( X )
%close all
%% LOAD MODEL AND DATA
Model = 'Macrophage_simple';
optModel = IQMmodel(strcat(Model,'.txt'));
IQMmakeMEXmodel(optModel,Model);
[paramNames] = IQMparameters(optModel);
[reactionNames] = IQMreactions(optModel);
initCon = IQMinitialconditions(optModel);
LoadData;
LoadDrugData;
TNF_Timepoints = DRUGDATA.timepoints;
TNF_Timepoints2 = DRUGDATA.timepoints2;
X = exp(X);
%% LOAD PARAMS
%values=load('allGoodValues.dat');
%ind=find(values(:,1)<min(values(:,1))*1.001,1);
%all_p = values(ind,:);
%X = all_p(2:end-3);
%scale_exp1 = all_p(end-2)
%scale_exp2 = all_p(end-1)
%scale_exp3 = all_p(end)
%% DEFINE COLORS
LPS0 = [0,0,0];
LPS100 = [0.5,0,0];
LPS1000 = [1,0,0];
dexa03 = [0,0,0.5];
dexa03_pre = [0,0,0.8];
dexa3 = [0,0.5,0];
dexa3_pre1 = [0,0.5,1];
dexa3_pre2 = [0.5,0.5,0.7];
dexa3_pre3 = [1,0.5,0.4];
%% SIMULATIONS
X = [X, 0, 0];
ss_simulation = IQMPsimulate(Model, 700, initCon, paramNames, X);
ss_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%%%%%%%%%% DEXAMETHASONE EXPERIMENTS FROM MARIA LIND MASTER THESIS %%%%%%
%ONLY LPS
X(end) = 0;
X(end-1) = 0;
LPS_0_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 10;
LPS_10_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 50;
LPS_50_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 100;
LPS_100_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 250;
LPS_250_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 500;
LPS_500_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 1000;
LPS_1000_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
figure(9)
errorbar([0,EXPDATA.lps], [EXPDATA.tnfaControl(end), EXPDATA.tnfa], [EXPDATA.semControl(end), EXPDATA.sem], 'ks','LineWidth', 1.5, 'MarkerFaceColor', 'k', 'MarkerSize', 4);
hold on
plot([0,EXPDATA.lps], [LPS_0_simulation.variablevalues(end,2),LPS_10_simulation.variablevalues(end,2),LPS_50_simulation.variablevalues(end,2),LPS_100_simulation.variablevalues(end,2),LPS_250_simulation.variablevalues(end,2),LPS_500_simulation.variablevalues(end,2),LPS_1000_simulation.variablevalues(end,2)], 'k', 'LineWidth', 2);
title('LPS dose-response');
xlabel('LPS (np/ml)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
figure(10)
hold on
errorbar(EXPDATA.timepoints, EXPDATA.tnfa1000, EXPDATA.sem1000,'s', 'Color', LPS1000, 'LineWidth', 1.5, 'MarkerFaceColor', LPS1000, 'MarkerSize', 4);
errorbar(EXPDATA.timepoints, EXPDATA.tnfa100, EXPDATA.sem100, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
errorbar(EXPDATA.timepoints, EXPDATA.tnfaControl, EXPDATA.semControl, 's', 'Color', LPS0, 'LineWidth', 1.5, 'MarkerFaceColor', LPS0, 'MarkerSize', 4);
plot(LPS_0_simulation.time, LPS_0_simulation.variablevalues(:,2), 'Color', LPS0, 'LineWidth', 2);
plot(LPS_100_simulation.time, LPS_100_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(LPS_1000_simulation.time, LPS_1000_simulation.variablevalues(:,2), 'Color', LPS1000, 'LineWidth', 2);
axis([0, 25, 0, 500]);
title('LPS');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('1000 uM LPS', '100 uM LPS', 'control', 'Location', 'northwest');
set(LEG,'FontSize',15);
%PRE-TREATMENT DEXA
X(end-1) = 0;
X(end) = 0.3;
Dexa_simulation_pre300 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, X);
X(end) = 3;
Dexa_simulation_pre3000 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, X);
%%addition of 100 ng/ml LPS
X(end-1) = 100;
X(end) = 0;
Dexa_simulation = IQMPsimulate(Model, 25, ss_simulation.statevalues(end,:), paramNames, X);
X(end) = 0.3;
Dexa_simulation_pre300.statevalues(end,end) = 0; %TNFa
Model_simulation2 = IQMPsimulate(Model, 25, Dexa_simulation_pre300.statevalues(end ,:), paramNames, X);
X(end) = 3;
Dexa_simulation_pre3000.statevalues(end,end) = 0; %TNFa
Model_simulation3 = IQMPsimulate(Model, 25, Dexa_simulation_pre3000.statevalues(end, :), paramNames, X);
X(end) = 0;
Model_simulation4 = IQMPsimulate(Model, 25, Dexa_simulation_pre300.statevalues(end,:), paramNames, X);
X(end) = 0;
Model_simulation5 = IQMPsimulate(Model, 25, Dexa_simulation_pre3000.statevalues(end, :), paramNames, X);
time_scale_exp2 = round(DRUGDATA.timepoints*1000/24)+1;
scale=lscov([Dexa_simulation.variablevalues(time_scale_exp2,2); Model_simulation2.variablevalues(time_scale_exp2,2); Model_simulation3.variablevalues(time_scale_exp2,2); Model_simulation4.variablevalues(time_scale_exp2,2); Model_simulation5.variablevalues(time_scale_exp2,2)], [DRUGDATA.tnfa100Control'; DRUGDATA.tnfaLowdex'; DRUGDATA.tnfaHighdex'; DRUGDATA.tnfaLowdex_pre'; DRUGDATA.tnfaHighdex_pre']);
figure(2);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfa100Control, DRUGDATA.sem100Control, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
hold on
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaLowdex_pre, DRUGDATA.semLowdex_pre, 's', 'Color', dexa03_pre, 'LineWidth', 1.5, 'MarkerFaceColor', dexa03_pre, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaHighdex_pre, DRUGDATA.semHighdex_pre, 's', 'Color', dexa3_pre3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaLowdex, DRUGDATA.semLowdex, 's', 'Color', dexa03, 'LineWidth', 1.5, 'MarkerFaceColor', dexa03, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaHighdex, DRUGDATA.semHighdex, 's', 'Color', dexa3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3, 'MarkerSize', 4);
plot(Dexa_simulation.time, scale.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(Model_simulation2.time, scale.*Model_simulation2.variablevalues(:,2),'Color', dexa03, 'LineWidth', 2);
plot(Model_simulation3.time, scale.*Model_simulation3.variablevalues(:,2),'Color', dexa3, 'LineWidth', 2);
plot(Model_simulation4.time, scale.*Model_simulation4.variablevalues(:,2),'Color', dexa03_pre, 'LineWidth', 2);
plot(Model_simulation5.time, scale.*Model_simulation5.variablevalues(:,2),'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 25, 0, 500]);
title('Dexamethasone + 100 uM LPS');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('control', '0.3 uM dexa (removed)', '3 uM dexa (removed)', '0.3 uM dexa', '3 uM dexa', 'Location', 'northwest');
set(LEG,'FontSize',15);
%%%%%%%%%% DEXAMETHASONE EXPERIMENTS FROM MARIA LIND MASTER THESIS: LONG-TIME EFFECT %%%%%%
X = X(1:(end-2));
X = [X, 0, 3]; %% Addition of 3 ug Dexa for 1 hour, no LPS
one_hour_sim = IQMPsimulate(Model, 1, ss_simulation.statevalues(end, :), paramNames, X);
one_hour_sim.statevalues(end,end) = 0; %TNFa = 0
X(end) = 0; %% 22 hours without dexa nor LPS
twentytwo_hour_sim = IQMPsimulate(Model, 22, one_hour_sim.statevalues(end,:), paramNames, X);
X(end-1) = 100; %% Trigger an inflammatory effect with LPS.
twentytwo_hour_sim.statevalues(end,end) = 0; %%TNFa = 0
pre1_sim = IQMPsimulate(Model, 7, twentytwo_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 0, 0];
sixteen_hour_sim = IQMPsimulate(Model, 16, one_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 100, 0];
sixteen_hour_sim.statevalues(end,end) = 0; %TNFa = 0
pre2_sim = IQMPsimulate(Model, 7, sixteen_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 100, 0];
pre3_sim = IQMPsimulate(Model, 7, one_hour_sim.statevalues(end,:), paramNames, X);
pre3_sim.variablevalues(round(TNF_Timepoints2 * 1000 / 7) + 1,2)
figure(3);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
hold on
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 4);
plot(Dexa_simulation.time, scale.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(pre1_sim.time, scale.*pre1_sim.variablevalues(:,2), 'Color', dexa3_pre1, 'LineWidth', 2);
plot(pre2_sim.time, scale.*pre2_sim.variablevalues(:,2), 'Color', dexa3_pre2, 'LineWidth', 2);
plot(pre3_sim.time, scale.*pre3_sim.variablevalues(:,2), 'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 7, 0, 210]);
title('Pre-treatments with dexa');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('control', 'pre1', 'pre2', 'pre3', 'Location', 'northwest');
set(LEG,'FontSize',15);
figure(89)
plot(1:7/1000:8, pre3_sim.statevalues(:,3), 'Color', dexa3_pre3, 'LineWidth', 2);
hold on
plot(16:7/1000:23, pre2_sim.statevalues(:,3), 'Color',dexa3_pre2, 'LineWidth', 2);
plot(22:7/1000:29, pre1_sim.statevalues(:,3), 'Color', dexa3_pre1, 'LineWidth', 2);
%plot(1:25/1000:26, Model_simulation5.statevalues(:,3), 'Color', dexa3_pre3, 'LineWidth', 2);
plot(0:0.001:1, one_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:25/1000:26, Model_simulation3.statevalues(:,3), 'Color', dexa3, 'LineWidth', 2);
figure(101)
statenames = IQMstates(Model);
for i = 1:9
subplot(3,3,i), plot(LPS_100_simulation.time, LPS_100_simulation.statevalues(:,i), 'Color', LPS100, 'LineWidth', 2);
hold on
subplot(3,3,i),plot(1:7/1000:8, pre3_sim.statevalues(:,i), 'Color', dexa3_pre3, 'LineWidth', 2);
subplot(3,3,i),plot(16:7/1000:23, pre2_sim.statevalues(:,i), 'Color',dexa3_pre2, 'LineWidth', 2);
subplot(3,3,i),plot(22:7/1000:29, pre1_sim.statevalues(:,i), 'Color', dexa3_pre1, 'LineWidth', 2);
%subplot(3,3,i),plot(1:25/1000:26, Model_simulation3.statevalues(:,i), 'Color', dexa3, 'LineWidth', 2);
subplot(3,3,i), plot(0:0.001:1, one_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(3,3,i),plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(3,3,i), plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
%subplot(3,3,i),plot(1:25/1000:26, Model_simulation2.statevalues(:,i), 'Color', dexa03, 'LineWidth', 2);
subplot(3,3,i),plot(1:25/1000:26, Model_simulation4.statevalues(:,i), 'Color', dexa03_pre, 'LineWidth', 2);
subplot(3,3,i),plot(1:25/1000:26, Model_simulation5.statevalues(:,i), 'Color', dexa3_pre3, 'LineWidth', 2);
title(statenames(i))
end
figure(102)
statenames = IQMstates(Model);
for i = 1:15
subplot(4,4,i), plot(LPS_100_simulation.time, LPS_100_simulation.reactionvalues(:,i), 'Color', LPS100, 'LineWidth', 2);
hold on
subplot(4,4,i),plot(1:7/1000:8, pre3_sim.reactionvalues(:,i), 'Color', dexa3_pre3, 'LineWidth', 2);
subplot(4,4,i),plot(16:7/1000:23, pre2_sim.reactionvalues(:,i), 'Color',dexa3_pre2, 'LineWidth', 2);
subplot(4,4,i),plot(22:7/1000:29, pre1_sim.reactionvalues(:,i), 'Color', dexa3_pre1, 'LineWidth', 2);
%subplot(4,4,i),plot(1:25/1000:26, Model_simulation3.reactionvalues(:,i), 'Color', dexa3, 'LineWidth', 2);
subplot(4,4,i), plot(0:0.001:1, one_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(4,4,i),plot(1:15/1000:16, sixteen_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(4,4,i), plot(1:21/1000:22, twentytwo_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
%subplot(4,4,i),plot(1:25/1000:26, Model_simulation2.reactionvalues(:,i), 'Color', dexa03, 'LineWidth', 2);
subplot(4,4,i),plot(1:25/1000:26, Model_simulation4.reactionvalues(:,i), 'Color', dexa03_pre, 'LineWidth', 2);
title(reactionNames(i))
end
a=1;
end
\ No newline at end of file
HypothesisB/LoadData.m
+ 4
− 13
View file @ 4497cdf2
......@@ -16,22 +16,21 @@ c4 = [13.15878 14.14623 18.20223]; l100_4 = [85.83585 89.46417 103.6841 105.3329
c6 = [14.30888 14.67227 18.2173]; l100_6 = [158.4887 154.0828 162.2422 255.6139 180.5526 190.212]; l1000_6 = [206.7315 213.3036 206.4169 348.8194 342.8387 306.1818];
c24 = [27.62139 34.96253 26.30299 57.61628 62.45974 69.22889]; l100_24 = [166.9165 147.1972 173.7321 630.5525 439.8181 446.5341]; l1000_24 = [212.4641 226.8367 232.186 558.759 568.0421 504.5358];
% Timepoints in the time response plot (figure 10)
% Timepoints in the time response plot
EXPDATA.timepoints = [0, 1, 2, 4, 6, 24];
% TNFA-response of the control in the time-response plot (figure 10)
% TNFA-response of the control in the time-response plot
EXPDATA.tnfaControl = [0, mean(c1), mean(c2), mean(c4), mean(c6), mean(c24)];
EXPDATA.semControl = [0, std(c1)/sqrt(3)/corrs3, std(c2)/sqrt(3)/corrs3, std(c4)/sqrt(3)/corrs3, std(c6)/sqrt(3)/corrs3, std(c24)/sqrt(6)/corrs6];
% TNFA-response of LPS conc. 100 ng/ml in the time-response plot (figure 10)
% TNFA-response of LPS conc. 100 ng/ml in the time-response plot
EXPDATA.tnfa100 = [0, mean(l100_1), mean(l100_2), mean(l100_4), mean(l100_6), mean(l100_24)];
EXPDATA.sem100 = [0, std(l100_1)/sqrt(3)/corrs3, std(l100_2)/sqrt(4)/corrs4, std(l100_4)/sqrt(6)/corrs6, std(l100_6)/sqrt(6)/corrs6, std(l100_24)/sqrt(6)/corrs6];
% TNFA-response of LPS conc. 1000 ng/ml in the time-response plot (figure 10)
% TNFA-response of LPS conc. 1000 ng/ml in the time-response plot
EXPDATA.tnfa1000 = [0, mean(l1000_1), mean(l1000_2), mean(l1000_4), mean(l1000_6), mean(l1000_24)];
EXPDATA.sem1000 = [0, std(l1000_1)/sqrt(3)/corrs3, std(l1000_2)/sqrt(6)/corrs6, std(l1000_4)/sqrt(6)/corrs6, std(l1000_6)/sqrt(6)/corrs6, std(l1000_24)/sqrt(6)/corrs6];
%x-axeln i figur 9 i Maria Linds master
EXPDATA.lps = [10,50,100,250,500,1000];
%LPS (ng/mL)
......@@ -42,18 +41,10 @@ c250 = [468.6055 382.417];
c500 = [465.6814 458.3753];
c1000 = [490.0401 485.2569];
%y-axeln i figur 9 i Maria Linds master
EXPDATA.tnfa = [mean(c10),mean(c50),mean(c100),mean(c250),mean(c500),mean(c1000)];
%residuals i figur 9 i Maria Linds master.
EXPDATA.sem = [std(c10)/sqrt(2)/corrs2,std(c50)/sqrt(2)/corrs2,std(c100)/sqrt(2)/corrs2,std(c250)/sqrt(2)/corrs2,std(c500)/sqrt(2)/corrs2,std(c1000)/sqrt(2)/corrs2];
%% Special SEM calculations
%EXPDATA.sem100 = max(EXPDATA.sem100./EXPDATA.tnfa100).*EXPDATA.tnfa100
%EXPDATA.sem1000 = max(EXPDATA.sem1000./EXPDATA.tnfa1000).*EXPDATA.tnfa1000
%EXPDATA.semControl = max(EXPDATA.semControl./EXPDATA.tnfaControl).*EXPDATA.tnfaControl
%EXPDATA.sem = max(EXPDATA.sem./EXPDATA.tnfa).*EXPDATA.tnfa
%% Special SEM calculations
EXPDATA.sem100 = [0,mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),mean(EXPDATA.sem100(2:end)),EXPDATA.sem100(end)]; %use SEM from dose-response (when higher)
EXPDATA.sem1000 = [0,mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),mean(EXPDATA.sem1000(2:end)),EXPDATA.sem1000(end)]; %even out SEM
......
HypothesisB/LoadDrugData.m
+ 16
− 48
View file @ 4497cdf2
%% DEXAMETHASONE DATA (INHIBITION OF TNF-alfa DATA). MARIA LINDH %%
%% DEXAMETHASONE DATA, MARIA LINDH %%
DRUGDATA = [];
% Timepoints in the drug response time-series (figure 12)
% Timepoints in the drug response time-series
DRUGDATA.timepoints = [0, 1, 2, 4, 6, 24];
%%% https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation
......@@ -22,25 +22,23 @@ c4 = [62.75599 57.82652]; l4 = [14.40909 17.0369 8.907761]; h4 = [2.258175 5.61
c6 = [194.3419 172.892]; l6 = [39.0532 43.15816 27.96148]; h6 = [16.71108 21.23619 12.19983]; l6p = [122.1991 197.0661 117.1118]; h6p = [50.12288 34.81426 32.16351];
c24 = [450.8582 485.3612]; l24 = [55.49924 53.36025 45.31882]; h24 = [29.17651 26.14517 24.43495]; l24p = [257.5461 391.8945 307.7493]; h24p = [96.68551 84.14325 97.12671];
% TNFA-response of the positive control (no drug) in the time-series (figure 12)
% TNFA-response of the positive control (no drug) in the time-series
DRUGDATA.tnfa100Control = [0, 0, mean(c2), mean(c4), mean(c6), mean(c24)];
DRUGDATA.sem100Control = [0, 0, std(c2)/sqrt(2)/corrs2, std(c4)/sqrt(2)/corrs2, std(c6)/sqrt(2)/corrs2, std(c24)/sqrt(2)/corrs2];
% TNFA-response of 0.3 �M dexamethasone in the time-series plot (figure 12)
% TNFA-response of 0.3 uM dexamethasone in the time-series plot
DRUGDATA.tnfaLowdex = [0, 0, 0, mean(l4), mean(l6), mean(l24)];
DRUGDATA.semLowdex = [0, 0, 0, std(l4)/sqrt(3)/corrs3, std(l6)/sqrt(3)/corrs3, std(l24)/sqrt(3)/corrs3];
% TNFA-response of 3 �M dexamethasone in the time-series plot (figure 12)
% TNFA-response of 3 uM dexamethasone in the time-series plot
DRUGDATA.tnfaHighdex = [0, 0, 0, mean(h4), mean(h6), mean(h24)];
DRUGDATA.semHighdex = [0, 0, 0, std(h4)/sqrt(2)/corrs2, std(h6)/sqrt(3)/corrs3, std(h24)/sqrt(3)/corrs3];
% TNFA-response of 0.3 �M dexamethasone (pre-treatment) in the time-series.
% (fig12)
% TNFA-response of 0.3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.tnfaLowdex_pre = [0, 0, 0, mean(l4p), mean(l6p), mean(l24p)];
DRUGDATA.semLowdex_pre = [0, 0, 0, std(l4p)/sqrt(3)/corrs3, std(l6p)/sqrt(3)/corrs3, std(l24p)/sqrt(3)/corrs3];
% TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig12)
% TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.tnfaHighdex_pre = [0, 0, 0, mean(h4p), mean(h6p), mean(h24p)];
DRUGDATA.semHighdex_pre = [0, 0, 0, std(h4p)/sqrt(3)/corrs3, std(h6p)/sqrt(3)/corrs3, std(h24p)/sqrt(3)/corrs3];
......@@ -51,7 +49,6 @@ DRUGDATA.semHighdex = [0,0,0,max(DRUGDATA.semHighdex),max(DRUGDATA.semHighdex),m
DRUGDATA.semLowdex_pre = [0,0,0,max(DRUGDATA.semLowdex_pre),max(DRUGDATA.semLowdex_pre),max(DRUGDATA.semLowdex_pre)];
DRUGDATA.semHighdex_pre = [0,0,0,max(DRUGDATA.semHighdex_pre),max(DRUGDATA.semHighdex_pre),max(DRUGDATA.semHighdex_pre)];
%% Fig 13 data (pre-treatment at different timepoints)
%Time of LPS incubation (hours) control control control pre-treatment 1 pre-treatment 1 pre-treatment 1 pre-treatment 2 pre-treatment 2 pre-treatment 2 pre-treatment 3 pre-treatment 3 pre-treatment 3
c0 = [2.47024 1.00613]; p1_0 = [0]; p2_0 = [0]; p3_0 = [0.656634];
c1 = [9.552152 9.736456 5.547666]; p1_1 = [2.686642]; p2_1 = [0]; p3_1 = [0];
......@@ -59,60 +56,31 @@ c2 = [27.86127 37.92371 23.21987]; p1_2 = [9.274983 8.05449 11.15372]; p2_2 = [4
c4 = [103.4877 172.0678 93.0803]; p1_4 = [65.6666 54.6981 55.74894]; p2_4 = [54.38446 51.23764 77.39728]; p3_4 = [17.42182 15.70416 35.74637];
c6 = [152.6309 220.9061 122.7121]; p1_6 = [89.44148 67.00409 85.37128]; p2_6 = [77.26416 86.26146 125.5514]; p3_6 = [32.67914 37.54742 61.43369];
%Annova
%Matr = [c6; p3_6; p2_6; p1_6]';
%[~,~,stats] = anova1(Matr);
%[c,~,~,gnames] = multcompare(stats,'CType','bonferroni');
% TIMEPOINTS
DRUGDATA.timepoints2 = [0, 1, 2, 4, 6];
% CONTROL: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% CONTROL: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre_control_tnfa = [mean(c0), mean(c1), mean(c2), mean(c4), mean(c6)];
DRUGDATA.pre_control_sem = [std(c0)/sqrt(2)/corrs2, std(c1)/sqrt(3)/corrs3, std(c2)/sqrt(3)/corrs3, std(c4)/sqrt(3)/corrs3, std(c6)/sqrt(3)/corrs3];
% PRETREATMENT1: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT1: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre1_tnfa = [0, 0, mean(p1_2), mean(p1_4), mean(p1_6)];
DRUGDATA.pre1_sem = [0, 0, std(p1_2)/sqrt(3)/corrs3, std(p1_4)/sqrt(3)/corrs3, std(p1_6)/sqrt(3)/corrs3];
% PRETREATMENT2: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT2: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre2_tnfa = [0, 0, mean(p2_2), mean(p2_4), mean(p2_6)];
DRUGDATA.pre2_sem = [0, 0, std(p2_2)/sqrt(3)/corrs3, std(p2_4)/sqrt(3)/corrs3, std(p2_6)/sqrt(3)/corrs3];
% PRETREATMENT3: TNFA-response of 3 �M dexamethasone (pre-treatment) in the time-series.
% (fig13)
% PRETREATMENT3: TNFA-response of 3 uM dexamethasone (pre-treatment) in the time-series.
DRUGDATA.pre3_tnfa = [0, 0, mean(p3_2), mean(p3_4), mean(p3_6)];
DRUGDATA.pre3_sem = [0, 0, std(p3_2)/sqrt(3)/corrs3, std(p3_4)/sqrt(3)/corrs3, std(p3_6)/sqrt(3)/corrs3];
%Set SEM to MAX(SEM)
DRUGDATA.pre_control_sem = [max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem),max(DRUGDATA.pre_control_sem)];
DRUGDATA.pre1_sem = [0,0,max(DRUGDATA.pre1_sem),max(DRUGDATA.pre1_sem),max(DRUGDATA.pre1_sem)];
DRUGDATA.pre2_sem = [0,0,max(DRUGDATA.pre2_sem),max(DRUGDATA.pre2_sem),max(DRUGDATA.pre2_sem)];
DRUGDATA.pre3_sem = [0,0,max(DRUGDATA.pre3_sem),max(DRUGDATA.pre3_sem),max(DRUGDATA.pre3_sem)];
%% Special SEM calculations
%sem1 = nanmean([DRUGDATA.sem100Control./DRUGDATA.tnfa100Control, DRUGDATA.semLowdex./DRUGDATA.tnfaLowdex, DRUGDATA.semHighdex./DRUGDATA.tnfaHighdex, DRUGDATA.semLowdex_pre./DRUGDATA.tnfaLowdex_pre, DRUGDATA.semHighdex_pre./DRUGDATA.tnfaHighdex_pre])
%DRUGDATA.sem100Control = sem1.*DRUGDATA.tnfa100Control;
%DRUGDATA.semLowdex = sem1.*DRUGDATA.tnfaLowdex;
%DRUGDATA.semHighdex = sem1.*DRUGDATA.tnfaHighdex;
%DRUGDATA.semLowdex_pre = sem1.*DRUGDATA.tnfaLowdex_pre;
%DRUGDATA.semHighdex_pre = sem1.*DRUGDATA.tnfaHighdex_pre;
%DRUGDATA.sem100Control = max(DRUGDATA.sem100Control./DRUGDATA.tnfa100Control).*DRUGDATA.tnfa100Control;
%DRUGDATA.semLowdex = max(DRUGDATA.semLowdex./DRUGDATA.tnfaLowdex).*DRUGDATA.tnfaLowdex;
%DRUGDATA.semHighdex = max(DRUGDATA.semHighdex./DRUGDATA.tnfaHighdex).*DRUGDATA.tnfaHighdex;
%DRUGDATA.semLowdex_pre = max(DRUGDATA.semLowdex_pre./DRUGDATA.tnfaLowdex_pre).*DRUGDATA.tnfaLowdex_pre;
%DRUGDATA.semHighdex_pre = max(DRUGDATA.semHighdex_pre./DRUGDATA.tnfaHighdex_pre).*DRUGDATA.tnfaHighdex_pre;
%sem2 = nanmean([DRUGDATA.pre_control_sem./DRUGDATA.pre_control_tnfa, DRUGDATA.pre1_sem./DRUGDATA.pre1_tnfa, DRUGDATA.pre2_sem./DRUGDATA.pre2_tnfa, DRUGDATA.pre3_sem./DRUGDATA.pre3_tnfa])
%DRUGDATA.pre_control_sem = sem2.*DRUGDATA.pre_control_tnfa;
%DRUGDATA.pre1_sem = sem2.*DRUGDATA.pre1_tnfa;
%DRUGDATA.pre2_sem = sem2.*DRUGDATA.pre2_tnfa;
%DRUGDATA.pre3_sem = sem2.*DRUGDATA.pre3_tnfa;
%DRUGDATA.pre_control_sem = max(DRUGDATA.pre_control_sem./DRUGDATA.pre_control_tnfa).*DRUGDATA.pre_control_tnfa;
%DRUGDATA.pre1_sem = max(DRUGDATA.pre1_sem./DRUGDATA.pre1_tnfa).*DRUGDATA.pre1_tnfa;
%DRUGDATA.pre2_sem = max(DRUGDATA.pre2_sem./DRUGDATA.pre2_tnfa).*DRUGDATA.pre2_tnfa;
%DRUGDATA.pre3_sem = max(DRUGDATA.pre3_sem./DRUGDATA.pre3_tnfa).*DRUGDATA.pre3_tnfa;
HypothesisB/Macrophage_simple.txt
+ 1
− 1
View file @ 4497cdf2
********** MODEL NAME
Macrophage_simple
Hypothesis B
********** MODEL NOTES
time in hours, concentrations in uM
......
HypothesisB/Macrophage_treatment.txt
+ 1
− 1
View file @ 4497cdf2
********** MODEL NAME
Macrophage_treatment
Treatment responses
********** MODEL NOTES
time in hours, concentrations in uM
......
HypothesisB/Optimization_simple.m
+ 12
− 38
View file @ 4497cdf2
clear all
close all
global EXPDATA
......@@ -10,7 +8,6 @@ global TNF_Timepoints
global DRUGDATA
global TNF_Timepoints2
global FID
global prev_cost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD THE MODEL
......@@ -27,51 +24,28 @@ TNF_Timepoints = DRUGDATA.timepoints;
TNF_Timepoints2 = DRUGDATA.timepoints2;
FID = fopen('allGoodValues.dat','wt');
prev_cost = 600;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THE OPTIMIZATION
psOpts = optimoptions(@particleswarm,'Display','iter');
saOpts = optimoptions(@simulannealbnd, 'HybridFcn',@fmincon,'Display','iter');%'TolFun', 1e-6, 'StallIterLimit', 1000,
%X = log(paramValues(1:14))';
X = [-2.5351 5.862 1.4977 -0.30429 -0.80806 -2.4145 5.2739 -4.1412 6.1595 6.9077 -6.9077 -6.8967 6.3553 6.6012 -1.9773 2.8764];
X = [-3.1396 6.1295 1.068 -0.94636 -0.84978 -2.4231 5.5636 -4.061 6.0711 7.0288 -7.1889 -7.0912 6.422 6.6677 -2.1251 2.8762];
nParams = 16;
%%
for i=1:20
i
%startGuess = X';
%startCost = cost_simple(startGuess');
%Define upper and lower bounds
lb = ones(16,1) * 1e-3;
%lb(5) = 0.5*60; %kon 0.5 uM/min = 30 uM/h
%lb(5) = 0.1; %kon !!!!
%lb(6) = 0.0001*60; %koff 0.0001 /min = 0.006 /h
%lb(6) = 0.0034*60*60; 12
lb=log(lb);
ub = ones(16,1) * 1e3;
%ub(5) = 1*60; %kon 60
%ub(6) = 0.01 * 60; %koff 0.6
%ub(6) = 0.007*60*60; 25
ub=log(ub);
format long
format compact
%Define upper and lower bounds
lb = ones(16,1) * 1e-3;
lb=log(lb);
ub = ones(16,1) * 1e3;
ub=log(ub);
format long
format compact
[optParamPS, minfunPS]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
[X, FVAL]=simulannealbnd(@cost_simple, optParamPS, lb,ub,saOpts);
%[X, FVAL]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
%[X, FVAL]=simulannealbnd(@cost_simple, X, lb,ub,saOpts);
save(sprintf('opt(%.2f), %s.mat',FVAL, datestr(now,'yymmdd-HHMMSS')),'X')
end
i
[optParamPS, minfunPS]=particleswarm(@cost_simple, nParams, lb,ub,psOpts);
[X, FVAL]=simulannealbnd(@cost_simple, optParamPS, lb,ub,saOpts);
save(sprintf('opt(%.2f), %s.mat',FVAL, datestr(now,'yymmdd-HHMMSS')),'X')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fclose(FID)
HypothesisB/cost_simple.m
+ 16
− 46
View file @ 4497cdf2
......@@ -6,24 +6,15 @@ function [cost] = cost_simple(param, shouldIPlot)
global Model
global TNF_Timepoints
global TNF_Timepoints2
global model1_EXPDATA
global DRUGDATA
global FID
global prev_cost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SIMULATION
param = exp(param);
%if param(6)/param(5) > 0.02 || param(6)/param(5) < 0.001
% cost = 500;
% return
%else
%param(end) = 0;
%%STEADY STATE SIMULATION (LPS = Dexa = 0)
param = [param, 0, 0];
try
......@@ -32,7 +23,7 @@ function [cost] = cost_simple(param, shouldIPlot)
ss_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%% LPS simulations from Maria Linds master thesis (figure 10 in her report).
%% LPS simulations
param(end) = 0;
param(end-1) = 0;
......@@ -52,38 +43,38 @@ LPS_1000_simulation = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.stateval
param(end-1) = 0;
doseresponseLPS = [LPS_10_simulation.variablevalues(end,2),LPS_50_simulation.variablevalues(end,2),LPS_100_simulation.variablevalues(end,2),LPS_250_simulation.variablevalues(end,2),LPS_500_simulation.variablevalues(end,2),LPS_1000_simulation.variablevalues(end,2)];
%% Dexamethasone simulations using data from Maria Linds master thesis (figure 12 in her report).
%% Dexamethasone simulations
% PRE-TREATMENT SIMULATIONS - macrophages are pre-treated with 0.3 or 3 �g
% PRE-TREATMENT SIMULATIONS - macrophages are pre-treated with 0.3 or 3 ug
% Dexa for 1 hour
param(end) = 0.3;
Dexa_simulation_pre300 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
param(end) = 3;
Dexa_simulation_pre3000 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
% 0 �M DEXAMETHASONE SIMULATION (positive control), only 100 ng/ml LPS
% 0 uM DEXAMETHASONE SIMULATION (positive control), only 100 ng/ml LPS
param = param(1:(end-2));
param = [param, 100, 0];
Dexa_simulation1 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end,:), paramNames, param);
% 0.3 �M DEXAMETHASONE SIMULATION
% 0.3 uM DEXAMETHASONE SIMULATION
param(end) = 0.3;
Dexa_simulation2 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end, :), paramNames, param);
% 3 �M DEXAMETHASONE SIMULATION
% 3 uM DEXAMETHASONE SIMULATION
param(end) = 3;
Dexa_simulation3 = IQMPsimulate(Model, TNF_Timepoints, ss_simulation.statevalues(end, :), paramNames, param);
% 0.3 �M DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0; % 1 % of 0.3 uM remains
% 0.3 uM DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0;
Dexa_simulation4 = IQMPsimulate(Model, TNF_Timepoints, Dexa_simulation_pre300.statevalues(end, :), paramNames, param);
% 3 �M DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0; % 1 % of 3 uM remains
% 3 uM DEXAMETHASONE SIMULATION (pre-treatment)
param(end) = 0;
Dexa_simulation5 = IQMPsimulate(Model, TNF_Timepoints, Dexa_simulation_pre3000.statevalues(end, :), paramNames, param);
%% Dexamethasone simulations using data from Maria Linds master thesis (figure 13 in her report).
%% Dexamethasone simulations
%Experiments with pre-treatment at different timepoints.
% % Control
......@@ -91,14 +82,14 @@ param = param(1:(end-2));
param = [param, 100, 0];
Pre_control_simulation = IQMPsimulate(Model, TNF_Timepoints2, ss_simulation.statevalues(end,:), paramNames, param);
% % 1 hour stimulation with 3 �g dexa (pre-treatment)
% % 1 hour stimulation with 3 ug dexa (pre-treatment)
param = param(1:(end-2));
param = [param, 0, 3];
one_hour_simulation = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, param);
one_hour_simulation.statevalues(end, end) = 0; %TNFa 0, buffer change
% % Pre-treatment 1
param(end) = 0; % 1 % of 3 uM remains
param(end) = 0;
twentytwo_hour_simulation = IQMPsimulate(Model, 22, one_hour_simulation.statevalues(end,:), paramNames, param);
param(end-1) = 100;
twentytwo_hour_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
......@@ -114,13 +105,6 @@ Pre2_simulation = IQMPsimulate(Model, TNF_Timepoints2, sixteen_hour_simulation.s
% % Pre-treatment 3
Pre3_simulation = IQMPsimulate(Model, TNF_Timepoints2, one_hour_simulation.statevalues(end,:), paramNames, param);
% % Pre-treatment long
%param(end-1) = 0;
%many_hour_simulation = IQMPsimulate(Model, 168, one_hour_simulation.statevalues(end,:), paramNames, param);
%param(end-1) = 100;
%many_hour_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%PreLong_simulation = IQMPsimulate(Model, TNF_Timepoints2, many_hour_simulation.statevalues(end,:), paramNames, param);
catch myError
cost = 1e6;
return
......@@ -138,21 +122,15 @@ Drug_simvalues4 = Dexa_simulation4.variablevalues(:,2);
Drug_simvalues5 = Dexa_simulation5.variablevalues(:,2);
%scale day-to-day differences between sets of data linearly
%scale_2=lscov([Control_simval; Pre1_simval; Pre2_simval; Pre3_simval], [DRUGDATA.pre_control_tnfa'; DRUGDATA.pre1_tnfa'; DRUGDATA.pre2_tnfa'; DRUGDATA.pre3_tnfa']);
scale=lscov([Drug_simvalues1; Drug_simvalues2; Drug_simvalues3; Drug_simvalues4; Drug_simvalues5], [DRUGDATA.tnfa100Control'; DRUGDATA.tnfaLowdex'; DRUGDATA.tnfaHighdex'; DRUGDATA.tnfaLowdex_pre'; DRUGDATA.tnfaHighdex_pre']);
%scale_exp3=lscov([LPS_0_simulation.variablevalues(:,2); LPS_100_simulation.variablevalues(:,2); LPS_1000_simulation.variablevalues(:,2)], [EXPDATA.tnfaControl'; EXPDATA.tnfa100'; EXPDATA.tnfa1000']);
%validation data
cost1 = sum((DRUGDATA.pre_control_tnfa - scale.*Control_simval').^2 ./ (DRUGDATA.pre_control_sem).^2);
cost2 = sum((DRUGDATA.pre1_tnfa(3:end) - scale.*Pre1_simval(3:end)').^2 ./ (DRUGDATA.pre1_sem(3:end)).^2);
cost3 = sum((DRUGDATA.pre2_tnfa(3:end) - scale.*Pre2_simval(3:end)').^2 ./ (DRUGDATA.pre2_sem(3:end)).^2);
cost4 = sum((DRUGDATA.pre3_tnfa(3:end) - scale.*Pre3_simval(3:end)').^2 ./ (DRUGDATA.pre3_sem(3:end)).^2);
%if abs(DRUGDATA.pre_control_tnfa(end) - scale.*PreLong_simulation.variablevalues(end,2)) > 20
% costextra = 50;
%else
% costextra = 0;
%end
%estimation data
cost5 = sum((DRUGDATA.tnfa100Control(3:end) - scale.*Drug_simvalues1(3:end)').^2 ./ (DRUGDATA.sem100Control(3:end)).^2);
cost6 = sum((DRUGDATA.tnfaLowdex(4:end) - scale.*Drug_simvalues2(4:end)').^2 ./ (DRUGDATA.semLowdex(4:end)).^2);
cost7 = sum((DRUGDATA.tnfaHighdex(4:end) - scale.*Drug_simvalues3(4:end)').^2 ./ (DRUGDATA.semHighdex(4:end)).^2);
......@@ -165,25 +143,17 @@ costLPS1000 = sum((EXPDATA.tnfa1000(2:end) - LPS_1000_simulation.variablevalues(
costLPSdose = sum((EXPDATA.tnfa - doseresponseLPS).^2 ./ (EXPDATA.sem).^2);
%size([DRUGDATA.pre_control_tnfa,DRUGDATA.pre2_tnfa(3:end),DRUGDATA.pre3_tnfa(3:end),DRUGDATA.tnfa100Control(3:end),DRUGDATA.tnfaLowdex(4:end),DRUGDATA.tnfaHighdex(4:end),DRUGDATA.tnfaHighdex_pre(4:end),EXPDATA.tnfa,EXPDATA.tnfa1000(2:end),EXPDATA.tnfa100(2:end),EXPDATA.tnfaControl(2:end),DRUGDATA.tnfaLowdex_pre(4:end)],2)
%costTNF = cost1 + cost2 + cost3 + cost4 + cost5 + cost6 + cost7 + cost8 + cost9;
costTNF = cost5 + cost6 + cost7 + cost8 + cost9;
costLPS = costLPS0 + costLPS100 + costLPS1000;
cost = costTNF + costLPSdose + costLPS; %+ costextra;
%cost_val = cost + cost2 + cost3
%cost = cost1 + cost2 + cost3 + cost4;
%cost1,cost2,cost3,cost4,cost5,cost6,cost7,cost8,cost9,costLPS0,costLPS100,costLPS1000,costLPSdose
%% CHI-SQUARE TEST
if cost < 83
%if cost < prev_cost
format short
prev_cost = cost;
fprintf(FID,'%4.10f %10.10f ',[cost, param(1:end-2), scale]); fprintf(FID,'\n');
end
%end
end
HypothesisB/plot_many_params.m
+ 16
− 185
View file @ 4497cdf2
%close all
clear all
close all
%% LOAD MODEL AND DATA
Model = 'Macrophage_simple';
optModel = IQMmodel(strcat(Model,'.txt'));
......@@ -13,8 +13,6 @@ LoadDrugData;
TNF_Timepoints = DRUGDATA.timepoints;
TNF_Timepoints2 = DRUGDATA.timepoints2;
FID = fopen('selectedValues.dat','wt');
%% DEFINE COLORS
LPS0 = [0,0,0];
LPS100 = [0.5,0,0];
......@@ -35,30 +33,8 @@ dexa03 = [213,94,0]./256; %red
dexa3_pre3 = [86,180,233]./256; %sky blue
dexa3 = [0,114,178]./256; %blue
dexa3_pre2 = [0,158,115]./256; %green
%dexa3_pre1 = [60,218,175]./256; %light green
%dexa3_pre1 = [80,238,195]./256; %light green
dexa3_pre1 = [80,0,0]./256; %brown
%% REDEFINE COLORS
%LPS0 = [0,0,0]; %black
%LPS1000 = [213,94,0]./256; %red
%dexa03 = [86,180,233]./256; %sky blue
%LPS100 = [204,121,167]./256; %purple pink
%dexa03_pre = [0,114,178]./256; %blue
%dexa3 = [240,228,66]./256; %yellow
%dexa3_pre3 = [230,159,0]./256; %orange
%dexa3_pre2 = [0,158,115]./256; %green
%dexa3_pre1 = [60,218,175]./256; %light green
%dexa3_pre1 = [80,238,195]./256; %light green
%dexa3_pre1 = [120,0,0]./256; %brown
%thres = 69; %all data chi2inv(0.95,51)
%thres = 65; %validation pretreatments chi2inv(0.95,48)
%thres = 52; %validation pretreatments chi2inv(0.95,37)
%thres = 80;
%% LOAD PARAMS
values1=load('allGoodValues1.dat'); %all data
values2=load('allGoodValues2.dat'); %all data
......@@ -71,10 +47,14 @@ values8=load('allGoodValues8.dat'); %all data
values9=load('allGoodValues9.dat'); %all data
values10=load('allGoodValues10.dat'); %no exp D
%Plot parameters that fit all data
%values = unique([values1;values2;values3;values6;values7;values8;values9],'rows');thres = 69;
%Plot parameters that fit estimation data
values = unique([values4;values5;values10],'rows');thres = 52;
ind=find(values(:,1)<thres);
res = 1000;
res = 3000;
all_param = ones(round(size(ind,1)/res), size(values1,2));
size(all_param)
......@@ -82,10 +62,8 @@ size(all_param)
%%
l=1;
for i = 1:res:size(ind, 1)
all_p = values(ind(i),:);
X = all_p(2:end-1);
......@@ -159,43 +137,11 @@ for i = 1:res:size(ind, 1)
X = X(1:(end-2));
X = [X, 100, 0];
pre3_sim = IQMPsimulate(Model, 7, one_hour_sim.statevalues(end,:), paramNames, X);
if 1.07*pre1_sim.variablevalues(end,2) < Dexa_simulation.variablevalues(round(1001/25*7),2)
best_pre1 = scale.*pre1_sim.variablevalues(:,2);
best_pre2 = scale.*pre2_sim.variablevalues(:,2);
best_pre3 = scale.*pre3_sim.variablevalues(:,2);
best_control = scale.*Dexa_simulation.variablevalues(:,2);
best_10 = LPS_10_simulation.variablevalues(end,2);
best_50 = LPS_50_simulation.variablevalues(end,2);
best_100 = LPS_100_simulation.variablevalues(round(1001/25*24),2);
best_250 = LPS_250_simulation.variablevalues(end,2);
best_500 = LPS_500_simulation.variablevalues(end,2);
best_1000 = LPS_1000_simulation.variablevalues(round(1001/25*24),2);
best_0 = LPS_0_simulation.variablevalues(:,2);
best_100_x = LPS_100_simulation.variablevalues(:,2);
best_1000_x = LPS_1000_simulation.variablevalues(:,2);
best_Dexa = scale.*Dexa_simulation.variablevalues(:,2);
best_sim2 = scale.*Model_simulation2.variablevalues(:,2);
best_sim3 = scale.*Model_simulation3.variablevalues(:,2);
best_sim4 = scale.*Model_simulation4.variablevalues(:,2);
best_sim5 = scale.*Model_simulation5.variablevalues(:,2);
%GR
basal = 0;%one_hour_sim.statevalues(1,7);
maxi = one_hour_sim.statevalues(end,7);
best_GR0 = (one_hour_sim.statevalues(:,7)-basal)./maxi;
best_GRbrown = (LPS_100_simulation.statevalues(:,7)-basal)./maxi;
best_GRpink = (pre3_sim.statevalues(:,7)-basal)./maxi;
best_GRgreen = (Model_simulation3.statevalues(:,7)-basal)./maxi;
end
% Uncomment to simulate only sustained response
%if 1.8*pre1_sim.variablevalues(end,2) < Dexa_simulation.variablevalues(round(1001/25*7),2)
%fprintf(FID,'%4.10f %10.10f ',all_p); fprintf(FID,'\n');
%all_p(6),all_p(7),all_p(8)
if l == 1
l = 2;
......@@ -253,7 +199,7 @@ for i = 1:res:size(ind, 1)
p3n_max = (pre3_sim.statevalues(:,8)-basal)./maxi;
%receptor
basal = 0;%one_hour_sim.statevalues(1,7);
basal = 0;
maxi = one_hour_sim.statevalues(end,7);
GR0_max = (one_hour_sim.statevalues(:,7)-basal)./maxi;
GR0_min = (one_hour_sim.statevalues(:,7)-basal)./maxi;
......@@ -333,7 +279,7 @@ for i = 1:res:size(ind, 1)
end
%receptor
basal = 0;%one_hour_sim.statevalues(1,7);
basal = 0;
maxi = one_hour_sim.statevalues(end,7);
for j = 1:size(GR0_min, 1)
......@@ -352,39 +298,7 @@ for i = 1:res:size(ind, 1)
GRgreen_max(j) = max(GRgreen_max(j), (Model_simulation3.statevalues(j,7)-basal)./maxi);
GRgreen_min(j) = min(GRgreen_min(j), (Model_simulation3.statevalues(j,7)-basal)./maxi);
end
%figure(105)
%hold on
%plot(0:0.001:1, (one_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%plot(22:7/1000:29, (pre1_sim.statevalues(:,9)-basal)./maxi, 'Color', dexa3_pre1, 'LineWidth', 0.5);
%plot(16:7/1000:23, (pre2_sim.statevalues(:,9)-basal)./maxi, 'Color',dexa3_pre2, 'LineWidth', 0.5);
%plot(1:7/1000:8, (pre3_sim.statevalues(:,9)-basal)./maxi, 'Color', dexa3_pre3, 'LineWidth', 0.5);
%plot(1:15/1000:16, (sixteen_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%plot(1:21/1000:22, (twentytwo_hour_sim.statevalues(:,9)-basal)./maxi, 'Color', LPS0, 'LineWidth', 0.5);
%figure(101)
%for i = 1:7
% basal = one_hour_sim.statevalues(1,i+5);
% %subplot(4,2,i), plot(LPS_100_simulation.time, LPS_100_simulation.statevalues(:,i+5), 'Color', LPS100, 'LineWidth', 0.5);
% hold on
% subplot(4,2,i),plot(1:7/1000:8, pre3_sim.statevalues(:,i+5)-basal, 'Color', dexa3_pre3, 'LineWidth', 0.5);
% subplot(4,2,i),plot(16:7/1000:23, pre2_sim.statevalues(:,i+5)-basal, 'Color',dexa3_pre2, 'LineWidth', 0.5);
% subplot(4,2,i),plot(22:7/1000:29, pre1_sim.statevalues(:,i+5)-basal, 'Color', dexa3_pre1, 'LineWidth', 0.5);
% %subplot(4,2,i),plot(1:25/1000:26, Model_simulation3.statevalues(:,i+5), 'Color', dexa3, 'LineWidth', 0.5);
% subplot(4,2,i), plot(0:0.001:1, one_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
% subplot(4,2,i),plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
% subplot(4,2,i), plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,i+5)-basal, 'Color', LPS0, 'LineWidth', 0.5);
%end
end
% figure(3)
% hold on
% plot(Dexa_simulation.time, scale_exp1.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 0.5);
% plot(pre1_sim.time, scale_exp1.*pre1_sim.variablevalues(:,2), 'Color', dexa3_pre1, 'LineWidth', 0.5);
% plot(pre2_sim.time, scale_exp1.*pre2_sim.variablevalues(:,2), 'Color', dexa3_pre2, 'LineWidth', 0.5);
% plot(pre3_sim.time, scale_exp1.*pre3_sim.variablevalues(:,2), 'Color', dexa3_pre3, 'LineWidth', 0.5);
%end
end
......@@ -397,8 +311,6 @@ errorbar(EXPDATA.lps, EXPDATA.tnfa, EXPDATA.sem, 'ks','LineWidth', 2.5, 'MarkerF
h = fill([EXPDATA.lps,flip(EXPDATA.lps)], [LPS_10_min,LPS_50_min,LPS_100_min,LPS_250_min,LPS_500_min,LPS_1000_min,LPS_1000_max,LPS_500_max,LPS_250_max,LPS_100_max,LPS_50_max,LPS_10_max], 'k');
set(h,'facealpha',.6,'EdgeColor','none')
%plot(EXPDATA.lps, [best_10,best_50,best_100,best_250,best_500,best_1000], 'Color', LPS0, 'LineWidth', 2);
title('A) LPS dose-response');
xlabel('LPS (ng/ml)');
ylabel('TNF (pg/mg tissue)');
......@@ -418,10 +330,6 @@ set(h,'facealpha',.6,'EdgeColor','none')
h = fill([LPS_0_simulation.time,flip(LPS_0_simulation.time)], ([LPS_0_min;flip(LPS_0_max)]'), LPS0);
set(h,'facealpha',.6,'EdgeColor','none')
%plot(LPS_0_simulation.time, best_0, 'Color', LPS0, 'LineWidth', 2);
%plot(LPS_100_simulation.time, best_100_x, 'Color', LPS100, 'LineWidth', 2);
%plot(LPS_1000_simulation.time, best_1000_x, 'Color', LPS1000, 'LineWidth', 2);
axis([0, 25, 0, 640]);
title('B) LPS time-response');
xlabel('Time (hours)');
......@@ -448,12 +356,6 @@ set(h,'facealpha',.6,'EdgeColor','none')
h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([sim5_min;flip(sim5_max)]'), dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
%plot(Dexa_simulation.time, best_Dexa, 'Color', LPS100, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim2, 'Color', dexa03, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim3, 'Color', dexa3, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim4, 'Color', dexa03_pre, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_sim5, 'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 25, 0, 640]);
title('C) Dexamethasone + LPS');
xlabel('Time (hours)');
......@@ -472,15 +374,8 @@ h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([GRbrown_min;flip(G
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([-1:0.001:0,flip(-1:0.001:0)], ([GR0_min;flip(GR0_max)]'), dexa3);
set(h,'facealpha',.6,'EdgeColor','none')
%plot(-1:0.001:0, best_GR0, 'Color', dexa3, 'LineWidth', 2);
%plot(-1:0.001:0, best_GR0, 'Color', dexa3, 'LineWidth', 2);
plot(Dexa_simulation.time, best_GRbrown, 'Color', LPS100, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_GRgreen, 'Color', dexa3, 'LineWidth', 2);
%plot(Dexa_simulation.time, best_GRpink, 'Color', dexa3_pre3, 'LineWidth', 2);
plot([0 0], [0 3], 'Color', LPS0, 'LineWidth', 0.5);
%xlim([-1 5])
axis([-1, 5, 0, 2.4]);
title('D) Mechanism of sustained response');
xlabel('Time (hours)');
......@@ -490,82 +385,22 @@ LEG = legend('3 uM Dexa', '3 uM Dexa (removed)', 'control (only LPS)', 'Location
set(LEG,'FontSize',14);
figure(2)
subplot(1,2,1)
hold on
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre1_min;flip(pre1_max)]'), dexa3_pre1);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre2_min;flip(pre2_max)]'), dexa3_pre2);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([Dexa_simulation.time,flip(Dexa_simulation.time)], ([Dexa_min;flip(Dexa_max)]'), LPS100);
set(h,'facealpha',.6,'EdgeColor','none')
h = fill([pre1_sim.time,flip(pre1_sim.time)], ([pre3_min;flip(pre3_max)]'), dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
plot(pre1_sim.time, best_pre1, 'Color', dexa3_pre1, 'LineWidth', 2);
plot(pre1_sim.time, best_pre2, 'Color', dexa3_pre2, 'LineWidth', 2);
plot(pre1_sim.time, best_pre3, 'Color', dexa3_pre3, 'LineWidth', 2);
plot(Dexa_simulation.time, best_control, 'Color', LPS100, 'LineWidth', 2);
axis([0, 7, 0, 260]);
title('Model predictions: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
set(gca,'FontSize',14);
%LEG = legend('control (only LPS)', '22 h wash before LPS', '16 h wash before LPS', 'dexa removed before LPS', 'Location', 'northwest');
%set(LEG,'FontSize',14);
subplot(1,2,2)
hold on
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 2.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 8);
axis([0, 7, 0, 260]);
title('Experimental data: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
ylabel('TNF (pg/mg tissue)');
set(gca,'FontSize',14);
LEG = legend('control (only LPS)', '22 h wash before LPS', '16 h wash before LPS', 'dexa removed before LPS', 'Location', 'northwest');
set(LEG,'FontSize',14);
%figure(101)
%for i = 1:7
% subplot(4,2,i),title(statenames(i+5))
%end
%scale2=lscov(mean([Dexa_min(DRUGDATA.timepoints2*40 + 1),Dexa_max(DRUGDATA.timepoints2*40 + 1)],2),DRUGDATA.pre_control_tnfa')
%scale2=0.84447
scale2=1
figure(102)
subplot(2,1,1)
hold on
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 2.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2 + 1, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2 + 17, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 8);
errorbar(DRUGDATA.timepoints2 + 23, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 2.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 8);
h = fill([0:0.025:6, flip(0:0.025:6)], scale2.*[Dexa_min(1:241);flip(Dexa_max(1:241))]', LPS100);
h = fill([0:0.025:6, flip(0:0.025:6)], [Dexa_min(1:241);flip(Dexa_max(1:241))]', LPS100);
set(h,'facealpha',.6,'EdgeColor','none')
%h=fill([0:0.001:1, flip(0:0.001:1)],[one_min; flip(one_max)]', LPS0);
%set(h,'facealpha',.6,'EdgeColor','none')
h=fill([1:7/1000:8, flip(1:7/1000:8)],scale2.*[pre3_min; flip(pre3_max)]', dexa3_pre3);
h=fill([1:7/1000:8, flip(1:7/1000:8)],[pre3_min; flip(pre3_max)]', dexa3_pre3);
set(h,'facealpha',.6,'EdgeColor','none')
h=fill([17:7/1000:24, flip(17:7/1000:24)],scale2.*[pre2_min; flip(pre2_max)]', dexa3_pre2);
h=fill([17:7/1000:24, flip(17:7/1000:24)],[pre2_min; flip(pre2_max)]', dexa3_pre2);
set(h,'facealpha',.6,'EdgeColor','none')
%h=fill([1:22/1000:23, flip(1:22/1000:23)],[twentytwo_min; flip(twentytwo_max)]', LPS0);
%set(h,'facealpha',.6,'EdgeColor','none')
h=fill([23:7/1000:30, flip(23:7/1000:30)],scale2.*[pre1_min; flip(pre1_max)]', dexa3_pre1);
h=fill([23:7/1000:30, flip(23:7/1000:30)],[pre1_min; flip(pre1_max)]', dexa3_pre1);
set(h,'facealpha',.6,'EdgeColor','none')
axis([0, 30, 0, max([Dexa_max(241),pre1_max(end)])+1]);
% hold on
% plot(0:0.001:1, (one_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
% plot(22:7/1000:29, (pre1_sim.statevalues(:,9)-basal)/maxi, 'Color', dexa3_pre1, 'LineWidth', 0.5);
% plot(16:7/1000:23, (pre2_sim.statevalues(:,9)-basal)/maxi, 'Color',dexa3_pre2, 'LineWidth', 0.5);
% plot(1:7/1000:8, (pre3_sim.statevalues(:,9)-basal)/maxi, 'Color', dexa3_pre3, 'LineWidth', 0.5);
%plot(1:15/1000:16, (sixteen_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
% plot(1:21/1000:22, (twentytwo_hour_sim.statevalues(:,9)-basal)/maxi, 'Color', LPS0, 'LineWidth', 0.5);
title('Model predictions: Dexamethasone, pre-treatments + LPS');
xlabel('Time (hours)');
......@@ -589,15 +424,11 @@ LEG = legend('control (only LPS)', 'dexa removed before LPS', '16 h wash before
set(LEG,'FontSize',14);
orderp = [1,3,4,8,7,2,11,10,9,14,15,12,13,5,6,16,17];
figure(103)
figure(3)
semilogy(1:size(all_param,2)-1,all_param(:,1+orderp),'-o')
%title('Acceptable parameter values');
%xlabel('Name of parameter');
set(gca,'XTick',1:size(all_param,2)-1)
set(gca,'XTickLabel',[paramNames(orderp(1:end-1));'scale'])
set(gca,'XTickLabelRotation',45)
ylabel('Parameter values');
set(gca,'FontSize',14);
axis([0, size(all_param,2), 9e-4, 2e3]);
fclose(FID)
\ No newline at end of file
HypothesisB/plot_treatment.m
+ 1
− 1
View file @ 4497cdf2
close all
clear all
%% LOAD MODEL AND DATA
Model = 'Macrophage_treatment';
optModel = IQMmodel(strcat(Model,'.txt'));
......
HypothesisB/simple_plot.m deleted 100644 → 0
+ 0
− 224
View file @ 91429cd6
function [ a ] = simple_plot( X )
%close all
%% LOAD MODEL AND DATA
Model = 'Macrophage_simple';
optModel = IQMmodel(strcat(Model,'.txt'));
IQMmakeMEXmodel(optModel,Model);
[paramNames] = IQMparameters(optModel);
[reactionNames] = IQMreactions(optModel);
initCon = IQMinitialconditions(optModel);
LoadData;
LoadDrugData;
TNF_Timepoints = DRUGDATA.timepoints;
TNF_Timepoints2 = DRUGDATA.timepoints2;
X = exp(X);
%% LOAD PARAMS
%values=load('allGoodValues.dat');
%ind=find(values(:,1)<min(values(:,1))*1.001,1);
%all_p = values(ind,:);
%X = all_p(2:end-3);
%scale_exp1 = all_p(end-2)
%scale_exp2 = all_p(end-1)
%scale_exp3 = all_p(end)
%% DEFINE COLORS
LPS0 = [0,0,0];
LPS100 = [0.5,0,0];
LPS1000 = [1,0,0];
dexa03 = [0,0,0.5];
dexa03_pre = [0,0,0.8];
dexa3 = [0,0.5,0];
dexa3_pre1 = [0,0.5,1];
dexa3_pre2 = [0.5,0.5,0.7];
dexa3_pre3 = [1,0.5,0.4];
%% SIMULATIONS
X = [X, 0, 0];
ss_simulation = IQMPsimulate(Model, 700, initCon, paramNames, X);
ss_simulation.statevalues(end, end) = 0; %TNFa 0 at experiment start
%%%%%%%%%% DEXAMETHASONE EXPERIMENTS FROM MARIA LIND MASTER THESIS %%%%%%
%ONLY LPS
X(end) = 0;
X(end-1) = 0;
LPS_0_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 10;
LPS_10_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 50;
LPS_50_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 100;
LPS_100_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 250;
LPS_250_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 500;
LPS_500_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
X(end-1) = 1000;
LPS_1000_simulation = IQMPsimulate(Model, 24, ss_simulation.statevalues(end,:), paramNames, X);
figure(9)
errorbar(EXPDATA.lps, EXPDATA.tnfa, EXPDATA.sem, 'ks','LineWidth', 1.5, 'MarkerFaceColor', 'k', 'MarkerSize', 4);
hold on
plot(EXPDATA.lps, [LPS_10_simulation.variablevalues(end,2),LPS_50_simulation.variablevalues(end,2),LPS_100_simulation.variablevalues(end,2),LPS_250_simulation.variablevalues(end,2),LPS_500_simulation.variablevalues(end,2),LPS_1000_simulation.variablevalues(end,2)], 'k', 'LineWidth', 2);
title('LPS dose-response');
xlabel('LPS (np/ml)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
figure(10)
hold on
errorbar(EXPDATA.timepoints, EXPDATA.tnfa1000, EXPDATA.sem1000,'s', 'Color', LPS1000, 'LineWidth', 1.5, 'MarkerFaceColor', LPS1000, 'MarkerSize', 4);
errorbar(EXPDATA.timepoints, EXPDATA.tnfa100, EXPDATA.sem100, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
errorbar(EXPDATA.timepoints, EXPDATA.tnfaControl, EXPDATA.semControl, 's', 'Color', LPS0, 'LineWidth', 1.5, 'MarkerFaceColor', LPS0, 'MarkerSize', 4);
plot(LPS_0_simulation.time, LPS_0_simulation.variablevalues(:,2), 'Color', LPS0, 'LineWidth', 2);
plot(LPS_100_simulation.time, LPS_100_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(LPS_1000_simulation.time, LPS_1000_simulation.variablevalues(:,2), 'Color', LPS1000, 'LineWidth', 2);
axis([0, 25, 0, 500]);
title('LPS');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('1000 uM LPS', '100 uM LPS', 'control', 'Location', 'northwest');
set(LEG,'FontSize',15);
%PRE-TREATMENT DEXA
X(end-1) = 0;
X(end) = 0.3;
Dexa_simulation_pre300 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, X);
X(end) = 3;
Dexa_simulation_pre3000 = IQMPsimulate(Model, 1, ss_simulation.statevalues(end,:), paramNames, X);
%%addition of 100 ng/ml LPS
X(end-1) = 100;
X(end) = 0;
Dexa_simulation = IQMPsimulate(Model, 25, ss_simulation.statevalues(end,:), paramNames, X);
X(end) = 0.3;
Dexa_simulation_pre300.statevalues(end,end) = 0; %TNFa
Model_simulation2 = IQMPsimulate(Model, 25, Dexa_simulation_pre300.statevalues(end ,:), paramNames, X);
X(end) = 3;
Dexa_simulation_pre3000.statevalues(end,end) = 0; %TNFa
Model_simulation3 = IQMPsimulate(Model, 25, Dexa_simulation_pre3000.statevalues(end, :), paramNames, X);
X(end) = 0;
Model_simulation4 = IQMPsimulate(Model, 25, Dexa_simulation_pre300.statevalues(end,:), paramNames, X);
X(end) = 0;
Model_simulation5 = IQMPsimulate(Model, 25, Dexa_simulation_pre3000.statevalues(end, :), paramNames, X);
time_scale_exp2 = round(DRUGDATA.timepoints*1000/24)+1;
scale=lscov([Dexa_simulation.variablevalues(time_scale_exp2,2); Model_simulation2.variablevalues(time_scale_exp2,2); Model_simulation3.variablevalues(time_scale_exp2,2); Model_simulation4.variablevalues(time_scale_exp2,2); Model_simulation5.variablevalues(time_scale_exp2,2)], [DRUGDATA.tnfa100Control'; DRUGDATA.tnfaLowdex'; DRUGDATA.tnfaHighdex'; DRUGDATA.tnfaLowdex_pre'; DRUGDATA.tnfaHighdex_pre']);
figure(2);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfa100Control, DRUGDATA.sem100Control, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
hold on
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaLowdex_pre, DRUGDATA.semLowdex_pre, 's', 'Color', dexa03_pre, 'LineWidth', 1.5, 'MarkerFaceColor', dexa03_pre, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaHighdex_pre, DRUGDATA.semHighdex_pre, 's', 'Color', dexa3_pre3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaLowdex, DRUGDATA.semLowdex, 's', 'Color', dexa03, 'LineWidth', 1.5, 'MarkerFaceColor', dexa03, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints, DRUGDATA.tnfaHighdex, DRUGDATA.semHighdex, 's', 'Color', dexa3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3, 'MarkerSize', 4);
plot(Dexa_simulation.time, scale.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(Model_simulation2.time, scale.*Model_simulation2.variablevalues(:,2),'Color', dexa03, 'LineWidth', 2);
plot(Model_simulation3.time, scale.*Model_simulation3.variablevalues(:,2),'Color', dexa3, 'LineWidth', 2);
plot(Model_simulation4.time, scale.*Model_simulation4.variablevalues(:,2),'Color', dexa03_pre, 'LineWidth', 2);
plot(Model_simulation5.time, scale.*Model_simulation5.variablevalues(:,2),'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 25, 0, 500]);
title('Dexamethasone + 100 uM LPS');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('control', '0.3 uM dexa (removed)', '3 uM dexa (removed)', '0.3 uM dexa', '3 uM dexa', 'Location', 'northwest');
set(LEG,'FontSize',15);
%%%%%%%%%% DEXAMETHASONE EXPERIMENTS FROM MARIA LIND MASTER THESIS: LONG-TIME EFFECT %%%%%%
X = X(1:(end-2));
X = [X, 0, 3]; %% Addition of 3 ug Dexa for 1 hour, no LPS
one_hour_sim = IQMPsimulate(Model, 1, ss_simulation.statevalues(end, :), paramNames, X);
one_hour_sim.statevalues(end,end) = 0; %TNFa = 0
X(end) = 0; %% 22 hours without dexa nor LPS
twentytwo_hour_sim = IQMPsimulate(Model, 22, one_hour_sim.statevalues(end,:), paramNames, X);
X(end-1) = 100; %% Trigger an inflammatory effect with LPS.
twentytwo_hour_sim.statevalues(end,end) = 0; %%TNFa = 0
pre1_sim = IQMPsimulate(Model, 7, twentytwo_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 0, 0];
sixteen_hour_sim = IQMPsimulate(Model, 16, one_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 100, 0];
sixteen_hour_sim.statevalues(end,end) = 0; %TNFa = 0
pre2_sim = IQMPsimulate(Model, 7, sixteen_hour_sim.statevalues(end,:), paramNames, X);
X = X(1:(end-2));
X = [X, 100, 0];
pre3_sim = IQMPsimulate(Model, 7, one_hour_sim.statevalues(end,:), paramNames, X);
pre3_sim.variablevalues(round(TNF_Timepoints2 * 1000 / 7) + 1,2)
figure(3);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre_control_tnfa, DRUGDATA.pre_control_sem, 's', 'Color', LPS100, 'LineWidth', 1.5, 'MarkerFaceColor', LPS100, 'MarkerSize', 4);
hold on
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre1_tnfa, DRUGDATA.pre1_sem, 's', 'Color', dexa3_pre1, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre1, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre2_tnfa, DRUGDATA.pre2_sem, 's', 'Color', dexa3_pre2, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre2, 'MarkerSize', 4);
errorbar(DRUGDATA.timepoints2, DRUGDATA.pre3_tnfa, DRUGDATA.pre3_sem, 's', 'Color', dexa3_pre3, 'LineWidth', 1.5, 'MarkerFaceColor', dexa3_pre3, 'MarkerSize', 4);
plot(Dexa_simulation.time, scale.*Dexa_simulation.variablevalues(:,2), 'Color', LPS100, 'LineWidth', 2);
plot(pre1_sim.time, scale.*pre1_sim.variablevalues(:,2), 'Color', dexa3_pre1, 'LineWidth', 2);
plot(pre2_sim.time, scale.*pre2_sim.variablevalues(:,2), 'Color', dexa3_pre2, 'LineWidth', 2);
plot(pre3_sim.time, scale.*pre3_sim.variablevalues(:,2), 'Color', dexa3_pre3, 'LineWidth', 2);
axis([0, 7, 0, 210]);
title('Pre-treatments with dexa');
xlabel('Time (h)');
ylabel('TNF-\alpha');
set(gca,'FontSize',18);
LEG = legend('control', 'pre1', 'pre2', 'pre3', 'Location', 'northwest');
set(LEG,'FontSize',15);
figure(89)
plot(1:7/1000:8, pre3_sim.statevalues(:,3), 'Color', dexa3_pre3, 'LineWidth', 2);
hold on
plot(16:7/1000:23, pre2_sim.statevalues(:,3), 'Color',dexa3_pre2, 'LineWidth', 2);
plot(22:7/1000:29, pre1_sim.statevalues(:,3), 'Color', dexa3_pre1, 'LineWidth', 2);
%plot(1:25/1000:26, Model_simulation5.statevalues(:,3), 'Color', dexa3_pre3, 'LineWidth', 2);
plot(0:0.001:1, one_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,3), 'Color', LPS0, 'LineWidth', 2);
plot(1:25/1000:26, Model_simulation3.statevalues(:,3), 'Color', dexa3, 'LineWidth', 2);
figure(101)
statenames = IQMstates(Model);
for i = 1:9
subplot(3,3,i), plot(LPS_100_simulation.time, LPS_100_simulation.statevalues(:,i), 'Color', LPS100, 'LineWidth', 2);
hold on
subplot(3,3,i),plot(1:7/1000:8, pre3_sim.statevalues(:,i), 'Color', dexa3_pre3, 'LineWidth', 2);
subplot(3,3,i),plot(16:7/1000:23, pre2_sim.statevalues(:,i), 'Color',dexa3_pre2, 'LineWidth', 2);
subplot(3,3,i),plot(22:7/1000:29, pre1_sim.statevalues(:,i), 'Color', dexa3_pre1, 'LineWidth', 2);
subplot(3,3,i),plot(1:25/1000:26, Model_simulation3.statevalues(:,i), 'Color', dexa3, 'LineWidth', 2);
subplot(3,3,i), plot(0:0.001:1, one_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(3,3,i),plot(1:15/1000:16, sixteen_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(3,3,i), plot(1:21/1000:22, twentytwo_hour_sim.statevalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(3,3,i),plot(1:25/1000:26, Model_simulation2.statevalues(:,i), 'Color', dexa03, 'LineWidth', 2);
subplot(3,3,i),plot(1:25/1000:26, Model_simulation4.statevalues(:,i), 'Color', dexa03_pre, 'LineWidth', 2);
title(statenames(i))
end
figure(102)
statenames = IQMstates(Model);
for i = 1:14
subplot(4,4,i), plot(LPS_100_simulation.time, LPS_100_simulation.reactionvalues(:,i), 'Color', LPS100, 'LineWidth', 2);
hold on
subplot(4,4,i),plot(1:7/1000:8, pre3_sim.reactionvalues(:,i), 'Color', dexa3_pre3, 'LineWidth', 2);
subplot(4,4,i),plot(16:7/1000:23, pre2_sim.reactionvalues(:,i), 'Color',dexa3_pre2, 'LineWidth', 2);
subplot(4,4,i),plot(22:7/1000:29, pre1_sim.reactionvalues(:,i), 'Color', dexa3_pre1, 'LineWidth', 2);
subplot(4,4,i),plot(1:25/1000:26, Model_simulation3.reactionvalues(:,i), 'Color', dexa3, 'LineWidth', 2);
subplot(4,4,i), plot(0:0.001:1, one_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(4,4,i),plot(1:15/1000:16, sixteen_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(4,4,i), plot(1:21/1000:22, twentytwo_hour_sim.reactionvalues(:,i), 'Color', LPS0, 'LineWidth', 2);
subplot(4,4,i),plot(1:25/1000:26, Model_simulation2.reactionvalues(:,i), 'Color', dexa03, 'LineWidth', 2);
subplot(4,4,i),plot(1:25/1000:26, Model_simulation4.reactionvalues(:,i), 'Color', dexa03_pre, 'LineWidth', 2);
title(reactionNames(i))
end
a=1;
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment