From 4497cdf2ade92877ae4a52e90ed0c0290b31c905 Mon Sep 17 00:00:00 2001
From: Elin Nyman <elin.nyman@liu.se>
Date: Wed, 5 Feb 2020 15:24:15 +0100
Subject: [PATCH] Update code
---
HypothesisA/LoadData.m | 15 +-
HypothesisA/LoadDrugData.m | 37 ++---
HypothesisA/Macrophage_simple.txt | 2 +-
HypothesisA/Optimization_simple.m | 60 ++-----
HypothesisA/cost_simple.m | 85 ++++------
HypothesisA/plot_many_params.m | 176 ++-------------------
HypothesisA/simple_plot.m | 226 ---------------------------
HypothesisB/LoadData.m | 17 +-
HypothesisB/LoadDrugData.m | 64 ++------
HypothesisB/Macrophage_simple.txt | 2 +-
HypothesisB/Macrophage_treatment.txt | 2 +-
HypothesisB/Optimization_simple.m | 50 ++----
HypothesisB/cost_simple.m | 62 ++------
HypothesisB/plot_many_params.m | 201 ++----------------------
HypothesisB/plot_treatment.m | 2 +-
HypothesisB/simple_plot.m | 224 --------------------------
16 files changed, 142 insertions(+), 1083 deletions(-)
delete mode 100644 HypothesisA/simple_plot.m
delete mode 100644 HypothesisB/simple_plot.m
diff --git a/HypothesisA/LoadData.m b/HypothesisA/LoadData.m
index 6f9379f..0d10fce 100644
--- a/HypothesisA/LoadData.m
+++ b/HypothesisA/LoadData.m
@@ -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
diff --git a/HypothesisA/LoadDrugData.m b/HypothesisA/LoadDrugData.m
index 0e1573b..056c313 100644
--- a/HypothesisA/LoadDrugData.m
+++ b/HypothesisA/LoadDrugData.m
@@ -1,8 +1,8 @@
-%% 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)];
diff --git a/HypothesisA/Macrophage_simple.txt b/HypothesisA/Macrophage_simple.txt
index f71d11e..eedfcc1 100644
--- a/HypothesisA/Macrophage_simple.txt
+++ b/HypothesisA/Macrophage_simple.txt
@@ -1,5 +1,5 @@
********** MODEL NAME
-Macrophage_simple
+Hypothesis A
********** MODEL NOTES
time in hours, concentrations in uM
diff --git a/HypothesisA/Optimization_simple.m b/HypothesisA/Optimization_simple.m
index 5fdee0b..a72b217 100644
--- a/HypothesisA/Optimization_simple.m
+++ b/HypothesisA/Optimization_simple.m
@@ -1,5 +1,4 @@
-
-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)
diff --git a/HypothesisA/cost_simple.m b/HypothesisA/cost_simple.m
index 6416812..c48d96c 100644
--- a/HypothesisA/cost_simple.m
+++ b/HypothesisA/cost_simple.m
@@ -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
diff --git a/HypothesisA/plot_many_params.m b/HypothesisA/plot_many_params.m
index 273e4a5..49016bc 100644
--- a/HypothesisA/plot_many_params.m
+++ b/HypothesisA/plot_many_params.m
@@ -1,5 +1,5 @@
-%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)
diff --git a/HypothesisA/simple_plot.m b/HypothesisA/simple_plot.m
deleted file mode 100644
index 605f7c4..0000000
--- a/HypothesisA/simple_plot.m
+++ /dev/null
@@ -1,226 +0,0 @@
-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
diff --git a/HypothesisB/LoadData.m b/HypothesisB/LoadData.m
index d344fd5..0d10fce 100644
--- a/HypothesisB/LoadData.m
+++ b/HypothesisB/LoadData.m
@@ -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
diff --git a/HypothesisB/LoadDrugData.m b/HypothesisB/LoadDrugData.m
index ac1b3a9..056c313 100644
--- a/HypothesisB/LoadDrugData.m
+++ b/HypothesisB/LoadDrugData.m
@@ -1,8 +1,8 @@
-%% 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;
diff --git a/HypothesisB/Macrophage_simple.txt b/HypothesisB/Macrophage_simple.txt
index 78d3bac..6b83ea8 100644
--- a/HypothesisB/Macrophage_simple.txt
+++ b/HypothesisB/Macrophage_simple.txt
@@ -1,5 +1,5 @@
********** MODEL NAME
-Macrophage_simple
+Hypothesis B
********** MODEL NOTES
time in hours, concentrations in uM
diff --git a/HypothesisB/Macrophage_treatment.txt b/HypothesisB/Macrophage_treatment.txt
index e7679d1..537d250 100644
--- a/HypothesisB/Macrophage_treatment.txt
+++ b/HypothesisB/Macrophage_treatment.txt
@@ -1,5 +1,5 @@
********** MODEL NAME
-Macrophage_treatment
+Treatment responses
********** MODEL NOTES
time in hours, concentrations in uM
diff --git a/HypothesisB/Optimization_simple.m b/HypothesisB/Optimization_simple.m
index de9cde2..dd0cf1f 100644
--- a/HypothesisB/Optimization_simple.m
+++ b/HypothesisB/Optimization_simple.m
@@ -1,5 +1,3 @@
-
-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)
diff --git a/HypothesisB/cost_simple.m b/HypothesisB/cost_simple.m
index 0117f5e..b1bdd29 100644
--- a/HypothesisB/cost_simple.m
+++ b/HypothesisB/cost_simple.m
@@ -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
diff --git a/HypothesisB/plot_many_params.m b/HypothesisB/plot_many_params.m
index 1be2d2b..4bb54a9 100644
--- a/HypothesisB/plot_many_params.m
+++ b/HypothesisB/plot_many_params.m
@@ -1,5 +1,5 @@
-%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
diff --git a/HypothesisB/plot_treatment.m b/HypothesisB/plot_treatment.m
index 59c9b3b..fb50a7c 100644
--- a/HypothesisB/plot_treatment.m
+++ b/HypothesisB/plot_treatment.m
@@ -1,5 +1,5 @@
close all
-clear all
+
%% LOAD MODEL AND DATA
Model = 'Macrophage_treatment';
optModel = IQMmodel(strcat(Model,'.txt'));
diff --git a/HypothesisB/simple_plot.m b/HypothesisB/simple_plot.m
deleted file mode 100644
index 7ad7e27..0000000
--- a/HypothesisB/simple_plot.m
+++ /dev/null
@@ -1,224 +0,0 @@
-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
--
GitLab