diff --git a/Readme.txt b/Readme.txt
index 3461013f2cd619d0e346bdf4cba685de0c692200..0b6b53ab569b2049bef063302e89e636f6df6469 100644
--- a/Readme.txt
+++ b/Readme.txt
@@ -1,4 +1,4 @@
-These scripts generate the results published in An organ-based multi-level model for glucose homeostasis: 
+These scripts generate the results published in An updated organ-based multi-level model for glucose homeostasis: 
 organ distributions, timing, and impact of blood flow.
 
 
@@ -9,16 +9,21 @@ Folders:
 ./data - contains all data files
 ./results - contains all acceptable parameters found
 ./scripts - contains all scripts except plotFigures
+./suport - contains neccessary toolboxes for model simulation
 
 Scripts:
-plotFigures.m - plots all figures in the manuscript except fig 8 (see 10.1074/jbc.M110.188987 for code to that figure)
+plotFigures.m - plots all figures in the manuscript.
 Setup.m - installs IQMtoolbox and compiles all model files. Used by plotFigures and runParallel.sh
-get_maxmin.m - get bounderies for estimation of parameter uncertainty.
-sample_params.m - samples all valid parameters found.
+runParallel.sh -top level script for running parameter estimation on the NSC super computer.
+main.m - script for running several parameter estimations in parallel.
+EstimateParameters.m EstimateParametersInputFat.m - script that runs parameter estimation, for the entire model and the adipocyte module respectively.
+simulateModel.m and simulateModelDalla.m - simulates glucose uptake variables and Dalla Man variables respectively.
+simulateSteadyState.m and simulateSteadyState.m - simulates the steady state of the models.
+costFun.m, costFunD.m, costFunFat.m, costFunVal.m, and costFunDVal.m- calculates the cost when finding optimal solution, (for the new glucose uptake variables or the Dalla Man variables), or when finding boundary of predictions.
+getMaxmin.m and getMaxminPred.m - get bounderies for estimation of parameter uncertainty for the fit to glucose uptake and for the prediction respectively (that is for plasma glucose and insulin).
+sampleParams.m - samples all valid parameters found.
+simInput.m - script that was used to simulate the input curve G_t used for simulation of adipocyte module (to get a good paramer range for these parameters).
+
 
-costFun.m - calculates cost, agreement between simulation and experimental data
-EstimateParameters.m - tries to find the best model parameters (only for model M4) with respect to agreement with data. 
-simulateModel - simulates and returns variables used for parameter estimation
-simulateSteadyState.m - simulates for a long period of time so that all states reaches steady state
 
 
diff --git a/plotFigures.m b/plotFigures.m
index 08b88f5776ead5d8e234fe1a8a6a64f973ea1fd4..526f0badfb860c35caa27e292003cdefda154a9b 100644
--- a/plotFigures.m
+++ b/plotFigures.m
@@ -1,7 +1,9 @@
-% plots all figures in manuscript
+
+
+%% plots all figures in manuscript
 
 close all
-clear all
+clear 
 
 cd('scripts')
 
@@ -10,7 +12,7 @@ Setup()
 
 %% fig 7 - fit to Brännmark data
 plotBrannmark 
-clear all
+clear 
 
 %% add necessary paths
 addpath('./Models')
@@ -22,12 +24,8 @@ addpath(genpath('./suport'))
 %% load data and parameters
 load('Coppack.mat');
 load('iozzo.mat');
-load('EGP.mat');
-load('Ra.mat');
-load('INS.mat');
-load('GLUT4m');
 
-load('./results/opt(12.92335).mat'); % optParam
+load('./results/optParam(24.549).mat'); % optParam
 load('./results/allAcceptableParams.mat'); % allParams
 
 %% nice colors
@@ -151,7 +149,7 @@ box on
 saveas(gcf, './results/fig2D.png')
 
 %% fig 3 - glucose uptake in adipose and muscle tissue, comparison model M2a and M2b
-model_names={'M1','M4m'};
+model_names={'M1','M4'};
 
 model=str2func(model_names{1});
 [~,paramsModA]=IQMparameters(model);
@@ -162,17 +160,15 @@ sim_muscleA = sim.variablevalues(:,ismember(IQMvariables(model),'U_idm'))';
 
 model=str2func(model_names{2}); 
 [pNames,paramsModB]=IQMparameters(model); 
-paramsFixed = [paramsModB(1:find(contains(pNames,'EGP1'))-1)' EGP Ra INS GLUT4m];
+paramsFixed = [paramsModB(1:find(contains(pNames,'p_bf')))'];
 paramsAllB = [paramsFixed optParam];
-[initsModB] = simulateSteadyState(model, paramsAllB);
+[initsModB] = simulateSteadyStated(model, paramsAllB);
 sim = model(sim_time, initsModB, paramsAllB);
 sim_adiposeB = sim.variablevalues(:,ismember(IQMvariables(model),'U_idf'))';
 sim_muscleB = sim.variablevalues(:,ismember(IQMvariables(model),'U_idm'))'; 
 
-
-sim = model(DATA.time, initsModB, paramsAllB);
-sim_adiposeDATA = sim.variablevalues(:,ismember(IQMvariables(model),'U_idf'))';
-sim_muscleDATA = sim.variablevalues(:,ismember(IQMvariables(model),'U_idm'))'; 
+% cost
+[sim_muscleDATA,sim_adiposeDATA,~,~] = simulateModel(DATA.time, model, optParam, paramsFixed);
 
 limit = chi2inv(0.95,numel(DATA.Muscle)+numel(DATA.Adipose))
 
@@ -183,7 +179,7 @@ cost = costMuscle + costAdipose
 
 % find min and max
 [sampledParams] = sampleParams(100,allParams);
-boundries=getMaxMin(sim_time, model, squeeze(sampledParams),paramsFixed);
+boundries=getMaxMin(sim_time, model, squeeze(sampledParams),paramsFixed,DATA,limit);
 
 % adipose
 figure()
@@ -241,7 +237,7 @@ saveas(gcf, './results/fig3Bmuscle.png')
 
 %% fig 4 - blood flow and glucose uptake adipose tissue
 
-model_name='M3m';
+model_name='M3';
 model=str2func(model_name);
 
 % normalize data
@@ -250,7 +246,7 @@ iozzo.clearance=100*iozzo.clearance/iozzo.clearance(1);
 
 % simulate experiments
 [pNames,paramsMod]=IQMparameters(model);
-paramsFixed = [paramsMod(1:find(contains(pNames,'EGP1'))-1)' EGP Ra INS GLUT4m];
+paramsFixed = [paramsMod(1:find(contains(pNames,'p_bf')))'];
 paramsAll=[paramsFixed optParam];
 [paramsI,paramsB,params]=deal(paramsAll);
 
@@ -260,24 +256,24 @@ paramsB(ismember(pNames,'p_bf')) = 5;
 paramsB(ismember(pNames,'INS_b')) = 0.8549; 
 paramsB(ismember(pNames,'bf_b')) = 3; 
 paramsBI = paramsB;
-paramsBI(ismember(pNames,'bradykinin')) = 10*paramsBI(ismember(pNames,'bradykinin')); % adding bradykinin, increases blood foow two-folds
+paramsBI(ismember(pNames,'bradykinin')) = 2.2*paramsBI(ismember(pNames,'bradykinin')); % adding bradykinin, increases blood foow two-folds
 
-[initsI] = simulateSteadyState(model, paramsI);
+[initsI] = simulateSteadyStated(model, paramsI);
 simI = model(sim_time, initsI, paramsI); 
 sim_uptakeI = simI.variablevalues(:,ismember(IQMvariables(model),'U_idf'))'; 
 AUCsim_uI = trapz(sim_time,sim_uptakeI);
 
-[inits] = simulateSteadyState(model, params);
+[inits] = simulateSteadyStated(model, params);
 sim = model(sim_time, inits, params); 
 sim_clearance = (sim.variablevalues(:,ismember(IQMvariables(model),'U_idf'))./sim.variablevalues(:,ismember(IQMvariables(model),'G')))'; 
 AUCsim_c = trapz(sim_time,sim_clearance);
 
-[initsB] = simulateSteadyState(model, paramsB);
+[initsB] = simulateSteadyStated(model, paramsB);
 simB = model(sim_time, initsB, paramsB); 
 sim_clearance_B = (simB.variablevalues(:,ismember(IQMvariables(model),'U_idf'))./sim.variablevalues(:,ismember(IQMvariables(model),'G')))'; 
 AUCsim_cB =trapz(sim_time,sim_clearance_B);
 
-[initsBI] = simulateSteadyState(model, paramsBI);
+[initsBI] = simulateSteadyStated(model, paramsBI);
 simBI = model(sim_time, initsBI, paramsBI); 
 sim_uptake_BI = simBI.variablevalues(:,ismember(IQMvariables(model),'U_idf'))';
 AUCsim_uB = trapz(sim_time,sim_uptake_BI);
@@ -344,12 +340,12 @@ saveas(gcf, './results/fig4.png')
 %% fig 5 - time-curves and total glucose uptake compared with data for all tissues
 
 clear model_name
-model_name='M4m';
+model_name='M4';
 model=str2func(model_name);
 
 [pNames,paramsMod]=IQMparameters(model); 
 pNames = pNames';
-paramsFixed = [paramsMod(1:find(contains(pNames,'EGP1'))-1)' EGP Ra INS GLUT4m];
+paramsFixed = [paramsMod(1:find(contains(pNames,'p_bf')))'];
 
 [sim_muscle,sim_adipose,sim_liver,sim_ii] = simulateModel(sim_time, model, optParam, paramsFixed);
 sim_tot = sim_muscle + sim_adipose + sim_liver + sim_ii;
@@ -358,9 +354,8 @@ sim_tot = sim_muscle + sim_adipose + sim_liver + sim_ii;
 figure()
 set(0,'DefaultLineLineWidth',3)
 set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 1]); % increase size of figure
-yticks([0:1:5])
+yticks([0:1:3])
 set(gca,'FontSize',20);
-axis([0 2*pi -1.5 5])
 
 % muscle
 ax1 = subplot(2,5,1);
@@ -368,9 +363,9 @@ plotArea(colorM,sim_muscle,sim_time); hold on
 plotSim(colorM,sim_muscle,sim_time); 
 plotErrorbar(colorM,DATA.time,DATA.Muscle,DATA.SEM1); 
 title('muscle tissue','FontSize',12);
-yticks([0:1:5])
+yticks([0:1:3])
 set(gca,'FontSize',20);
-axis([0 400 -0.5 5])
+axis([0 400 -0.5 3.5])
 
 % adipose
 ax2 = subplot(2,5,2);
@@ -378,28 +373,28 @@ plotArea(colorA,sim_adipose,sim_time); hold on
 plotSim(colorA,sim_adipose,sim_time); 
 plotErrorbar(colorA,DATA.time,DATA.Adipose,DATA.SEM2); 
 title('adipose tissue','FontSize',20);
-yticks([0:1:5])
+yticks([0:1:3])
 set(gca,'FontSize',20);
-axis([0 400 -0.5 5])
+axis([0 400 -0.5 3.5])
 
 % liver
 ax3 = subplot(2,5,3);
 plotArea(colorL,sim_liver,sim_time); hold on
 plotSim(colorL,sim_liver,sim_time); 
 title('liver','FontSize',20);
-yticks([0:1:5])
+yticks([0:1:3])
 set(gca,'FontSize',20);
-axis([0 400 -0.5 5])
+axis([0 400 -0.5 3.5])
 
 % insulin independent uptake 
 ax4 = subplot(2,5,4);
 plotArea(colorII,sim_ii,sim_time); hold on 
 plotSim(colorII,sim_ii,sim_time); 
 title('other','FontSize',20);
-yticks([0:1:5])
+yticks([0:1:3])
 set(gca,'FontSize',20);
 
-axis([0 400 -0.5 5])
+axis([0 400 -0.5 3.5])
 
 % total uptake 
 load('GU_high.mat');
@@ -431,40 +426,87 @@ ylabel('mg/kg');
 set(gca,'xtick',[])
 set(gca,'FontSize',20);
 title('Total postprandial uptake','FontSize',20);
-legend('framed: simulation AUC [0-360 min]','filled: data AUC [0-360 min]','Location', 'northwest')
+legend('framed: simulation AUC [0-360 min]','filled: data total uptake (Gerich 2000)','Location', 'northwest')
 set(gca,'Xtick',1:13,'XTickLabel',{'             muscle',[],[],'               adipose',[],[],'              liver',[],[],'               brain + kidneys',[],[],'total',[],[]});
 
 saveas(gcf, './results/fig5.png')
 
-% C - Validation, plasma insulin and glucose
+%% fig 5C - Validation, plasma insulin and glucose
+
+global SEM_EGP
+global EGP_dalla
+global SEM_S
+global S_dalla
+global SEM_I
+global pi_dalla
+global SEM_G
+global pg_dalla
+global SEM_Ra
+global Ra_dalla
+global SEM_U
+global GU_dalla
+
+load('EGP_dalla.mat');
+load('EGP_high.mat');
+load('EGP_low.mat');
+
+load('S_dalla.mat');
+load('S_high.mat');
+load('S_low.mat');
+
+load('pi_dalla.mat');
+load('pi_high.mat');
+load('pi_low.mat');
+
+load('pg_dalla.mat');
+load('pg_high.mat');
+load('pg_low.mat');
+
+load('Ra_dalla.mat');
+load('Ra_high.mat');
+load('Ra_low.mat');
+
+load('GU_dalla.mat');
+load('GU_high.mat');
+load('GU_low.mat');
+
+load('SEM_EGP.mat');
+load('SEM_S.mat');
+load('SEM_I.mat');
+load('SEM_G.mat');
+load('SEM_Ra.mat');
+load('SEM_U.mat');
+
 
 load Adams.mat OGTT
+model_name='M4pred';
+model=str2func(model_name);
+[pNames,paramsMod]=IQMparameters(model); 
+pNames = pNames';
+paramsFixed = [paramsMod(1:find(contains(pNames,'p_bf')))'];
+OGTT.Glucose = OGTT.Glucose(1:end-1); % removing the last datapoint glucose
+OGTT.SEMg = OGTT.SEMg(1:end-1);
+OGTT.Insulin = OGTT.Insulin; 
 [~,I] = mink(OGTT.SEMg(2:end),1); 
 OGTT.SEMg(I+1) = mean(OGTT.SEMg(2:end)); % changing the smallest SEM to that of the mean SEM
+[~,I] = mink(OGTT.SEMi(2:end),1); 
+OGTT.SEMi(I+1) = mean(OGTT.SEMi(2:end)); % changing the smallest SEM to that of the mean SEM
 
-model_name='M1';
-model=str2func(model_name);
+limit_95 = chi2inv(0.95,numel(OGTT.Glucose(2:end))+numel(OGTT.Insulin(2:end))); % The first data point is removed since it is used as initial value
 
-[pNames,paramsMod]=IQMparameters(model); 
-[inits] = simulateSteadyStated(model, paramsMod);
-inits(28) = OGTT.Glucose(1)*1.88; % G_p
-inits(31) = OGTT.Insulin(1)*0.05; % I_p
 sim_time = [OGTT.time(1):0.1:OGTT.time(end)];
-sim = model(OGTT.time, inits, paramsMod'); 
-IpVal = sim.variablevalues(:,ismember(IQMvariables(model),'I'))'; 
-GpVal = sim.variablevalues(:,ismember(IQMvariables(model),'G'))'; 
-
-simPlot = model(sim_time, inits, paramsMod');  
-Gp = simPlot.variablevalues(:,ismember(IQMvariables(model),'G'))'; 
-Ip = simPlot.variablevalues(:,ismember(IQMvariables(model),'I'))'; 
-
-% cost
-limit_95 = chi2inv(0.95,numel(OGTT.Glucose(2:end-1))+numel(OGTT.Insulin(2:end))) % The first data point is removed since it is used as initial value
-
-cost1 = nansum((((OGTT.Glucose(2:end-1)-GpVal(2:end-1)').^2)./(OGTT.SEMg(2:end-1).^2)))
-cost2 = nansum((((OGTT.Insulin(2:end)-IpVal(2:end)').^2)./(OGTT.SEMi(2:end).^2)))
+if isfile('./results/optParamsPred.mat')
+    load 'boundriesVal.mat' boundriesVal
+    load 'optParamsPred.mat' optParamsPred
+else
+[boundriesVal,optParamsPred,costPred]=getMaxMinPred(sim_time, model, squeeze(allParams(:,1:16)),paramsFixed,OGTT);
+end
+paramsAll = [paramsFixed optParamsPred];
 
-cost = cost1 + cost2
+[inits] = simulateSteadyStated(model, paramsAll);
+sim = model(sim_time, inits, paramsAll); 
+Ip = sim.variablevalues(:,ismember(IQMvariables(model),'I'))'; 
+Gp = sim.variablevalues(:,ismember(IQMvariables(model),'G'))'; 
 
 % plot
 colorData = [0, 0, 0]; % total
@@ -475,70 +517,43 @@ set(0,'DefaultLineLineWidth',3)
 
 % plasma glucose
 subplot(1,2,1);
-plotErrorbar(colorData,OGTT.time(1:end-1),OGTT.Glucose(1:end-1),OGTT.SEMg(1:end-1)); hold on
+plotUncertainty(colorSim,boundriesVal,sim_time,'g'); hold on
+plotErrorbar(colorData,OGTT.time(1:end-1),OGTT.Glucose(1:end),OGTT.SEMg(1:end)); 
 plotSim(colorSim,Gp,sim_time); 
 title('Plasma glucose','FontSize',15);
+xlim([0 250])
 set(gca,'FontSize',15);
 ylabel('mg/dl');
-legend('data','simulation','Location', 'northeast')
+legend('uncertainty','data','simulation','Location', 'northeast')
 set(gca,'linewidth',1.4)
 
-
 % plasma insulin
 subplot(1,2,2);
+plotUncertainty(colorSim,boundriesVal,sim_time,'i'); hold on
 plotErrorbar(colorData,OGTT.time,OGTT.Insulin,OGTT.SEMi); hold on
 plotSim(colorSim,Ip,sim_time); 
 title('Plasma Insulin','FontSize',15);
-% yticks([0:1:3])
-% ylim([0 3.5])
+xlim([0 250])
 ylabel('pmol/L');
 set(gca,'FontSize',15);
 set(gca,'linewidth',1.4)
-legend('data','simulation','Location', 'northeast')
-
+legend('uncertainty','data','simulation','Location', 'northeast')
 
 saveas(gcf, './results/GlucoseInsulin.png')
 
-
 %% fig 6 - fit to Dalla Man data
-load('EGP_ori.mat');
-load('EGP_dalla.mat');
-load('EGP_high.mat');
-load('EGP_low.mat');
-
-load('S_ori.mat');
-load('S_dalla.mat');
-load('S_high.mat');
-load('S_low.mat');
-
-load('pi_ori.mat');
-load('pi_dalla.mat');
-load('pi_high.mat');
-load('pi_low.mat');
-
-load('pg_ori.mat');
-load('pg_dalla.mat');
-load('pg_high.mat');
-load('pg_low.mat');
-
-load('Ra_ori.mat');
-load('Ra_dalla.mat');
-load('Ra_high.mat');
-load('Ra_low.mat');
-
-load('GU_ori.mat');
-load('GU_dalla.mat');
-load('GU_high.mat');
-load('GU_low.mat');
 
 % simulate
-clear model_name
-model_name='M1';
+sim_time = [pg_dalla(1,1):0.01:pg_dalla(end,1)];
+model_name='M4';
 model=str2func(model_name);
+
 [pNames,paramsMod]=IQMparameters(model); 
-sim_time = [pg_dalla(1,1):0.01:pg_dalla(end,1)];
-[inits] = simulateSteadyStated(model, paramsMod);
-simModel = model(sim_time, inits, paramsMod); 
+pNames = pNames';
+paramsFixed = [paramsMod(1:find(contains(pNames,'p_bf')))'];
+paramsAll = [paramsFixed optParam];
+[inits] = simulateSteadyStated(model, paramsAll);
+simModel = model(sim_time, inits, paramsAll); 
 
 % plot
 figure()
@@ -548,7 +563,6 @@ subplot(3,2,1)
 plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'G')),'k-');
 hold on
 plot(pg_dalla(:,1),pg_dalla(:,2),'k.')
-hold on
 plotUncertaintyD(pg_high,pg_low,sim_time)
 ylabel('(mg/dl)')
 title('Plasma Glucose')
@@ -557,7 +571,6 @@ subplot(3,2,2)
 plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'I')),'k-');
 hold on
 plot(pi_dalla(:,1),pi_dalla(:,2),'k.')
-hold on
 plotUncertaintyD(pi_high,pi_low,sim_time)
 ylabel('(pmol/L)')
 title('Plasma Insulin')
@@ -566,7 +579,6 @@ subplot(3,2,3)
 plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'EGP')),'k-');
 hold on
 plot(EGP_dalla(:,1),EGP_dalla(:,2),'k.')
-hold on 
 plotUncertaintyD(EGP_high,EGP_low,sim_time)
 hold off
 ylabel('(mg/kg/min)')
@@ -576,7 +588,6 @@ subplot(3,2,4)
 plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'Ra')),'k-');
 hold on
 plot(Ra_dalla(:,1),Ra_dalla(:,2),'k.')
-hold on
 plotUncertaintyD(Ra_high,Ra_low,sim_time)
 ylabel('(mg/kg/min)')
 title('Glucose Rate of Appearance')
@@ -585,13 +596,10 @@ subplot(3,2,6)
 plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'S')),'k-');
 hold on
 plot(S_dalla(:,1),S_dalla(:,2),'k.')
-hold on
 plotUncertaintyD(S_high,S_low,sim_time)
 ylabel('(pmol/kg/min)')
 title('Insulin secretion S')
 
-model_name='M4m';
-model=str2func(model_name);
 [sim_muscle,sim_adipose,sim_liver,sim_ii] = simulateModel(sim_time, model, optParam, paramsFixed);
 sim_tot = sim_muscle + sim_adipose + sim_liver + sim_ii;
 
@@ -599,10 +607,10 @@ subplot(3,2,5)
 plot(simModel.time, sim_tot,'k-');
 hold on
 plot(GU_dalla(:,1),GU_dalla(:,2),'k.')
-hold on
 plotUncertaintyD(GU_high,GU_low,sim_time)
 title('Glucose utilization')
 
+%%
 cd ..
 %----------------functions-----------------------
 
diff --git a/scripts/EstimateParameters.m b/scripts/EstimateParameters.m
index e7fd8f0c20da1a2c2a8061dea0492ffd4b5d061b..4350cc99bdcd90b52c49519aa8fe2af6dfc49ca4 100644
--- a/scripts/EstimateParameters.m
+++ b/scripts/EstimateParameters.m
@@ -1,4 +1,18 @@
 function [] =  EstimateParameters(worker,now)
+
+global SEM_EGP
+global EGP_dalla
+global SEM_S
+global S_dalla
+global SEM_I
+global pi_dalla
+global SEM_G
+global pg_dalla
+global SEM_Ra
+global Ra_dalla
+global SEM_U
+global GU_dalla
+
 % estimates parameters, runs from file it is placed in
  s=rng('shuffle'); %do not note remove this line
  s=rng(s.Seed+worker); %do not note remove this line
@@ -8,27 +22,35 @@ addpath('./Models')
 addpath(genpath('./suport/'))
 addpath(genpath('./results'))
 
-% initial stuff
-load ('Coppack.mat'); 
-load('EGP.mat')
-load('Ra.mat')
-load('INS.mat')
-load('GLUT4m.mat')
+% load data
+load('Coppack.mat'); 
 
-model_name='M4m';
+load('SEM_EGP.mat');
+load('EGP_dalla.mat');
+load('SEM_S.mat');
+load('S_dalla.mat');
+load('SEM_I.mat');
+load('pi_dalla.mat');
+load('SEM_G.mat');
+load('pg_dalla.mat');
+load('SEM_Ra.mat');
+load('Ra_dalla.mat');
+load('SEM_U.mat');
+load('GU_dalla.mat');
 
 % make model
+model_name='M4';
 model=str2func(model_name);
 
 % parameters
 [pNames,paramsMod]=IQMparameters(model); 
-
-paramsFixed = [paramsMod(1:find(contains(pNames,'EGP1'))-1)' EGP' Ra' INS' GLUT4m'];
-startParam = paramsMod(length(paramsFixed)+1:end);
-nparams=length(startParam);
+load 'optParamF.mat' optParam; % optimal parameters from fat module
+paramsFixed = paramsMod(1:find(contains(pNames,'p_bf')))';
+startParam = [paramsMod(length(paramsFixed)+1:end-length(optParam))' optParam];
+nparams = length(startParam);
 optNames = pNames(length(paramsFixed)+1:end);
 
-limit=chi2inv(0.95,numel(DATA.Adipose)+numel(DATA.Muscle));
+limit = chi2inv(0.95,numel(DATA.Muscle)+numel(DATA.Adipose)+numel(pg_dalla(:,1))*6-11);
 folder=sprintf('./results/%s/%s/',model_name,now);
 
 if ~exist(folder,'dir')
@@ -39,12 +61,12 @@ end
 
 options_ps=optimoptions(@particleswarm, 'Display','iter');
 options_s=optimoptions(@simulannealbnd, 'Display','iter','HybridFcn',@fmincon);
-ub = log(startParam.*10000);
-lb = log(startParam./10000);
+ub = log(startParam.*10);
+lb = log(startParam./10);
 
 fid=fopen(sprintf('%s/validParams-%s-%i.csv', folder, model_name,s.Seed),'at');
 
-costfunc = @(p) costFun(DATA,model,p,fid,limit,paramsFixed);
+costfunc = @(p) costFunD(DATA,model,p,fid,limit,paramsFixed);
 
 optParam_ps = startParam;
 [optParam_ps, mincost_ps] = particleswarm(costfunc, nparams, lb, ub, options_ps);
diff --git a/scripts/EstimateParametersInput.m b/scripts/EstimateParametersInputFat.m
similarity index 64%
rename from scripts/EstimateParametersInput.m
rename to scripts/EstimateParametersInputFat.m
index 410566dd397aecd9ad0c32cae72508e38f0e6b6e..72e245f146b083bb4d5f94d539815b1c11158964 100644
--- a/scripts/EstimateParametersInput.m
+++ b/scripts/EstimateParametersInputFat.m
@@ -1,21 +1,20 @@
-function [] =  EstimateParametersInput(worker,now)
+function [] =  EstimateParametersInputFat(worker,now)
 % estimates parameters, runs from file it is placed in
- s=rng('shuffle'); %do not note remove this line
- s=rng(s.Seed+worker); %do not note remove this line
+ s=rng('shuffle'); 
+ s=rng(s.Seed+worker); 
 
 addpath('./data')
 addpath('./Models')
 addpath(genpath('./results'))
-addpath(genpath('./../MATLAB/IQMtools/'))
+addpath(genpath('./suport/IQMtools/'))
 
 % load data
 load ('Coppack.mat'); 
-load('EGP.mat')
-load('Ra.mat')
+
 load('INS.mat')
-load('GLUT4m.mat')
+load('G_t.mat')
 
-model_name={'M4m'}; 
+model_name={'M4fat'}; 
 
 for k=1:length(model_name)
 
@@ -24,12 +23,12 @@ model=str2func(model_name{k});
 % params
 [pNames,paramsMod]=IQMparameters(model); 
 
-paramsFixed = [paramsMod(1:find(contains(pNames,'EGP1'))-1)' EGP Ra INS GLUT4m]; % G_t'
+paramsFixed = [paramsMod(1:find(contains(pNames,'INS1'))-1)' INS G_t]; 
 startParam = paramsMod(length(paramsFixed)+1:end)';
 nparams=length(startParam);
 optNames = pNames(length(paramsFixed)+1:end);
 
-limit=chi2inv(0.95,numel(DATA.Adipose)+numel(DATA.Muscle));
+limit=chi2inv(0.95,numel(DATA.Adipose));
 folder=sprintf('./results/%s/%s/',model_name{k},now);
 
 if ~exist(folder,'dir')
@@ -41,23 +40,25 @@ end
 options_ps=optimoptions(@particleswarm, 'Display','iter');
 options_s=optimoptions(@simulannealbnd, 'Display','iter','HybridFcn',@fmincon);
 
-ub = log(startParam.*10000);
-lb = log(startParam./10000);
+ub = log(startParam.*1000);
+lb = log(startParam./1000);
+
 
 fid=fopen(sprintf('%s/validParams-%s-%i.csv', folder, model_name{k},s.Seed),'at');
 
-costfunc = @(p) costFun(DATA,model,p,fid,limit,paramsFixed);
+costfunc = @(p) costFunFat(DATA,model,p,fid,limit,paramsFixed);
 
 [optParam_ps, mincost_ps] = particleswarm(costfunc, nparams, lb, ub, options_ps);
 [optParam, mincost] = simulannealbnd(costfunc,optParam_ps, lb, ub, options_s);
 fclose(fid);
 
 % Save results
-save(sprintf('%s/opt(%.5f)-%i.mat',folder,mincost,s.Seed) ,'optParam')
+optParam = exp(optParam);
+save(sprintf('%s/optF(%.5f)-%i.mat',folder,mincost,s.Seed) ,'optParam')
 
 if numel(dir(sprintf('%s/optNames.mat',folder)))==0
     save(sprintf('%s/optNames.mat',folder),'optNames')
 end
     
 end
-end
\ No newline at end of file
+end
diff --git a/scripts/Models/M0.txt b/scripts/Models/M0.txt
index 69758184b00e78b9f949462cbd1b6d6730fb26b7..b0fb47944b94c829e6beaec5a77996223efcb6a1 100644
--- a/scripts/Models/M0.txt
+++ b/scripts/Models/M0.txt
@@ -1,9 +1,7 @@
 ********** MODEL NAME
 M0
-
 ********** MODEL NOTES
 
-
 ********** MODEL STATES
 d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
 d/dt(IRp)               = v1basal+v1c-v1d-v1g   %%2
@@ -13,25 +11,26 @@ d/dt(IRi)               = v1e-v1r               %%5
 d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
 d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
 d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
-d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
 d/dt(X)                 = v3b-v3a               %%10
-d/dt(Xp)                = v3a-v3b               %%12
-d/dt(PKB)               = -v4a+v4b+v4h          %%13   
-d/dt(PKB308p)           = v4a-v4b-v4c           %%14
-d/dt(PKB473p)           = -v4e+v4f-v4h          %%15
-d/dt(PKB308p473p)       = v4c+v4e-v4f           %%16
-d/dt(mTORC1)            = v5b-v5a               %%17         
-d/dt(mTORC1a)           = v5a-v5b               %%18
-d/dt(mTORC2)            = -v5c+v5d              %%19
-d/dt(mTORC2a)           = v5c-v5d               %%20
-d/dt(AS160)             = v6b1-v6f1             %%21
-d/dt(AS160p)            = v6f1-v6b1             %%22
-d/dt(GLUT4m)            = (v7f-v7b)             %%23
-d/dt(GLUT4)             = -v7f+v7b              %%24
-d/dt(S6K)               = v9b1-v9f1             %%26
-d/dt(S6Kp)              = v9f1-v9b1             %%27
-d/dt(S6)                = v9b2-v9f2             %%28
-d/dt(S6p)               = v9f2-v9b2             %%29                                               
+d/dt(Xp)                = v3a-v3b               %%11
+d/dt(PKB)               = -v4a+v4b+v4h          %%12   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%13
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%14
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%15
+d/dt(mTORC1)            = v5b-v5a               %%16         
+d/dt(mTORC1a)           = v5a-v5b               %%17
+d/dt(mTORC2)            = -v5c+v5d              %%18
+d/dt(mTORC2a)           = v5c-v5d               %%19
+d/dt(AS160)             = v6b1-v6f1             %%20
+d/dt(AS160p)            = v6f1-v6b1             %%21
+d/dt(GLUT4m)            = v7f-v7b               %%22
+d/dt(GLUT4)             = -v7f+v7b              %%23
+d/dt(S6K)               = v9b1-v9f1             %%24
+d/dt(S6Kp)              = v9f1-v9b1             %%25
+d/dt(S6)                = v9b2-v9f2             %%26
+d/dt(S6p)               = v9f2-v9b2             %%27
+                                              
 d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
 d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                          
 d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                        
@@ -179,25 +178,24 @@ k_e2 = 339  %%87
 D = 78000  %%88    
    
 ********** MODEL VARIABLES
+I = I_p/V_I                                                                      
+G = G_p/V_G  
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
-EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_mmax = (1-part)*(V_m0+V_mX*INS)                                                
-V_fmax = part*(V_f0+V_fX*INS)                                                 
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                 
 E = 0                                                                            
-S = gamma*I_po                                                                                                                       
-I = I_p/V_I                                                                      
-G = G_p/V_G                                                                      
+S = gamma*I_po                                                                                                                                                                                                          
 HE = (-m_5*S)+m_6                                                                
 m_3 = HE*m_1/(1-HE)                                                              
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
 k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)     
+V_mmax = (1-part)*(V_m0+V_mX*INS) 
 U_idm = V_mmax*G_t/(K_m0+G_t)
-U_idf = scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3)                                                
-U_id = U_idm+U_idf                                                               
+U_idf = scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3)                                                                                                            
 U = U_ii+U_id                                                                    
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
+U_id = U_idm+U_idf 
 
 ********** MODEL REACTIONS                                                                                                          
 v1a      = IR*k1a*(INS+5)*1e-3
@@ -238,9 +236,6 @@ v9b1     = S6Kp*k9b1
 v9f2     = S6*k9f2*S6Kp
 v9b2     = S6p*k9b2
 
-glucoseuptake = k8*GLUT4m*(G_t*5/170) + glut1*(G_t*5/170) + kbf*((INS+5)*1e-3)
-
-
 ********** MODEL FUNCTIONS
 
 
diff --git a/scripts/Models/M1.txt b/scripts/Models/M1.txt
index 5236f8094be2b1cac7ea59634a028698a5f013d7..f7788fee4c34dc66299c10e734c07b5395bb62b2 100644
--- a/scripts/Models/M1.txt
+++ b/scripts/Models/M1.txt
@@ -32,6 +32,7 @@ d/dt(S6K)               = v9b1-v9f1             %%26
 d/dt(S6Kp)              = v9f1-v9b1             %%27
 d/dt(S6)                = v9b2-v9f2             %%28
 d/dt(S6p)               = v9f2-v9b2             %%29                                               
+
 d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
 d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                        
 d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S                                         
@@ -73,7 +74,8 @@ GLUT4(0)            = 100
 S6K(0)              = 100
 S6Kp(0)             = 0
 S6(0)               = 100
-S6p(0)              = 0                               
+S6p(0)              = 0 
+                              
 G_p(0) = 178                                                                   
 G_t(0) = 130                                                                     
 I_l(0) = 4.5                                                                 
@@ -164,33 +166,27 @@ k_p3 = 0.009  %%70
 k_p4 = 0.0618  %%71                                                              
 k_i = 0.0079  %%72                                                               
 
-V_m0 = 2.5    %%74                                                                   
-V_mX = 0.047  %%75                                                                
-K_m0 = 225.59  %%76                                             
-V_f0 = 2.5  %%77                                                                 
-V_fX = 0.047  %%78                                                               
-K_f0 = 225.59  %%79                                                              
-p_2U = 0.0331  %%80                                                              
-part = 0.2  %%81                                                                 
-K = 2.3  %%82                                                                    
-alpha = 0.05  %%83                                                               
-beta = 0.11  %%84                                                                
-gamma = 0.5  %%85                                                                
-k_e1 = 0.0005  %%86                                                              
-k_e2 = 339  %%87                                                                 
-D = 78000  %%88    
-
-K_l0 = 225.59
-V_l0 = 2.5                                                                       
-V_lX = 0.047  
+V_m = 2.5                                                                       
+V_mX = 0.047                                                                  
+K_m = 225.59                                                               
+p_2U = 0.0331                                                                 
+K = 2.3                                                                     
+alpha = 0.05                                                                
+beta = 0.11                                                                 
+gamma = 0.5                                                                  
+k_e1 = 0.0005                                                               
+k_e2 = 339                                                                  
+D = 78000     
+part_m = 0.5 
+part_f = 0.44 
+part_l = 1.67 
+U_ii = 0.8447
 
 ********** MODEL VARIABLES
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
 EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_mmax = (1-part)*(V_m0+V_mX*INS)                                                
-V_fmax = part*(V_f0+V_fX*INS) 
-V_lmax = (1-part)*(V_l0+V_lX*G_t)                                                   
+V_mmax = V_m+V_mX*INS                                                     
 Ee = k_e1*(G_p - k_e2)                                                                            
 S = gamma*I_po                                                                                                                       
 I = I_p/V_I                                                                      
@@ -200,13 +196,12 @@ m_3 = HE*m_1/(1-HE)
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
 k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)     
-U_idm = 0.5*V_mmax*G_t/(K_m0+G_t)
-U_idf = 0.5*scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3)                                                                                                              
+U_idm = part_m*V_mmax*G_t/(K_m+G_t)
+U_idf = part_f*scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3) 
+U_idl = part_l*U_idm                                                                                                          
 U = U_ii+U_id                                                                    
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
-U_liver = 45*U_idm/27
-U_ii = 23*U_liver/45
-U_id = U_idm+U_idf+U_liver
+U_id = U_idm+U_idf+U_idl
 
 ********** MODEL REACTIONS                                                                                                          
 v1a      = IR*k1a*(INS+5)*1e-3
diff --git a/scripts/Models/M2a.txt b/scripts/Models/M2a.txt
index a640722cc7cf31a0eb50a2b2bb80e315a9d4bb2e..af1d728f02e8e79f7a6251c557fda9a438ae5fd4 100644
--- a/scripts/Models/M2a.txt
+++ b/scripts/Models/M2a.txt
@@ -2,7 +2,7 @@
 M2a
 
 ********** MODEL NOTES
-
+M2a - different parameters for interstitial insulin breakdown and transport
 
 ********** MODEL STATES
 d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
@@ -13,41 +13,40 @@ d/dt(IRi)               = v1e-v1r               %%5
 d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
 d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
 d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
-d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
 d/dt(X)                 = v3b-v3a               %%10
-d/dt(Xp)                = v3a-v3b               %%12
-d/dt(PKB)               = -v4a+v4b+v4h          %%13   
-d/dt(PKB308p)           = v4a-v4b-v4c           %%14
-d/dt(PKB473p)           = -v4e+v4f-v4h          %%15
-d/dt(PKB308p473p)       = v4c+v4e-v4f           %%16
-d/dt(mTORC1)            = v5b-v5a               %%17         
-d/dt(mTORC1a)           = v5a-v5b               %%18
-d/dt(mTORC2)            = -v5c+v5d              %%19
-d/dt(mTORC2a)           = v5c-v5d               %%20
-d/dt(AS160)             = v6b1-v6f1             %%21
-d/dt(AS160p)            = v6f1-v6b1             %%22
-d/dt(GLUT4m)            = (v7f-v7b)             %%23
-d/dt(GLUT4)             = -v7f+v7b              %%24
-d/dt(S6K)               = v9b1-v9f1             %%26
-d/dt(S6Kp)              = v9f1-v9b1             %%27
-d/dt(S6)                = v9b2-v9f2             %%28
-d/dt(S6p)               = v9f2-v9b2             %%29                                             
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                            
-d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                       
-d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S                                         
+d/dt(Xp)                = v3a-v3b               %%11
+d/dt(PKB)               = -v4a+v4b+v4h          %%12   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%13
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%14
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%15
+d/dt(mTORC1)            = v5b-v5a               %%16         
+d/dt(mTORC1a)           = v5a-v5b               %%17
+d/dt(mTORC2)            = -v5c+v5d              %%18
+d/dt(mTORC2a)           = v5c-v5d               %%19
+d/dt(AS160)             = v6b1-v6f1             %%20
+d/dt(AS160p)            = v6f1-v6b1             %%21
+d/dt(GLUT4m)            = v7f-v7b               %%22
+d/dt(GLUT4)             = -v7f+v7b              %%23
+d/dt(S6K)               = v9b1-v9f1             %%24
+d/dt(S6Kp)              = v9f1-v9b1             %%25
+d/dt(S6)                = v9b2-v9f2             %%26
+d/dt(S6p)               = v9f2-v9b2             %%27             
+       
+                                                      
+d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
+d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                      
+d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S  
+d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                            
 d/dt(Q_sto1) = -k_gri*Q_sto1                                                     
 d/dt(Q_sto2) = (-k_empt*Q_sto2)+k_gri*Q_sto1                                     
 d/dt(Q_gut) = (-k_abs*Q_gut)+k_empt*Q_sto2                                       
 d/dt(I_1) = -k_i*(I_1-I)                                                         
-d/dt(I_d) = -k_i*(I_d-I_1)                                                       
-d/dt(INS) = (-p_2U*INS)+p_2U*(I-I_b)                                             
+d/dt(I_d) = -k_i*(I_d-I_1)                                                           
 d/dt(I_po) = (-gamma*I_po)+S_po                                                  
-d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                
-d/dt(INS_new) = (-k1*INS_new)+k2*(I-I_b)
-d/dt(INS_f) = (-C1_f*INS_f)+C2_f*(I-I_b)
-d/dt(INS_f_bf) = (-C1_bf*INS_f_bf)+C2_bf*(I-I_b)
-
+d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                        
+d/dt(INS_f) = (-p_2U*INS_f)+p_2U*(I-I_b)  
+d/dt(INS) = V_1-V_2    
 
 IR(0)               = 100
 IRp(0)              = 0
@@ -75,7 +74,8 @@ GLUT4(0)            = 100
 S6K(0)              = 100
 S6Kp(0)             = 0
 S6(0)               = 100
-S6p(0)              = 0                               
+S6p(0)              = 0    
+             
 G_p(0) = 178                                                                   
 G_t(0) = 130                                                                     
 I_l(0) = 4.5                                                                 
@@ -85,12 +85,12 @@ Q_sto2(0) = 0
 Q_gut(0) = 0                                                                     
 I_1(0) = 25                                                                      
 I_d(0) = 25                                                                      
-INS(0) = 0                                                                     
 I_po(0) = 3.6                                                  
 Y(0) = 0
-INS_new(0) = 0
+INS_f(0) = 0
+INS(0) = 0
 
-********** MODEL PARAMETERS
+********** MODEL PARAMETERS 
 
 diabetes = 1  %%1                                                                
 k1a = 0.633141  %%2                                                              
@@ -137,8 +137,8 @@ k9b2 = 30.9967  %%42
 km9 = 5872.68  %%43                                                              
 n9 = 0.985466  %%44                                                              
 kbf = 1e+06  %%45                                                                
-scaleModel = 2.1e-06  %%46  
-                                                     
+scaleModel = 2.1e-06  %%46   
+
 V_G = 1.88  %%47                                                                 
 k_1 = 0.065  %%48                                                                
 k_2 = 0.079  %%49                                                                
@@ -159,19 +159,11 @@ k_gri = 0.0558  %%63
 f = 0.9  %%64                                                                    
 b = 0.82  %%65                                                                   
 d = 0.01  %%66                                                                   
-BW = 78  %%67                                                                    
 k_p1 = 2.7  %%68                                                                 
 k_p2 = 0.0021  %%69                                                              
 k_p3 = 0.009  %%70                                                               
 k_p4 = 0.0618  %%71                                                              
 k_i = 0.0079  %%72                                                               
-U_ii = 1.47  %%73                                                                
-V_l0 = 2.5     %% 74                                                                  
-V_mX = 0.047  %%75                                                               
-K_l0 = 225.59  %%76                                                              
-V_f0 = 2.5  %%77                                                                 
-V_fX = 0.047  %%78                                                               
-K_f0 = 225.59  %%79                                                              
 p_2U = 0.0331  %%80                                                              
 part = 0.2  %%81                                                                 
 K = 2.3  %%82                                                                    
@@ -179,72 +171,56 @@ alpha = 0.05  %%83
 beta = 0.11  %%84                                                                
 gamma = 0.5  %%85                                                                
 k_e1 = 0.0005  %%86                                                              
-k_e2 = 339  %%87                                                                 
-D = 78000  %%88    
-
-scale10 = 0.5 %% 89
-Ip_b = 7 %% 90
-be = 3 %% 91
-
-k2 = 2.15 %% 90
-k1 = 1.3 %% 91
-C1_f = 1.25 %% 92
-C2_f = 0.6 %% 93
-V_m0 = 2.5 %% 94
-V_mx_new = 0.047 %% 96
-K_m0 = 225.59 %% 98
-BRF = 1 %% 100
-bf_b = 3 %% 101
-INS_b = 0.8549 %% 102
-scale_bf = 0.05e2 %% 103
-C1_bf = 0.075 %% 104
-C2_bf = 0.06 %% 105
-scale1 = 1 %% 106
-scale2 = 0.13 %% 107
-scale3 = 2.4 %% 108
-scale4 = 0.4 %% 109
-scale6 = 1.5 %% 110
-scale7 = 0.7 %% 111
-scale8 = 1.8 %% 112
-scale9 = 14 %% 113
-k_gluin = 1300 %% 114
-k_G6P = 15 %% 115
-V_G6Pmax = 90 %% 116
-pl = 0.5 %% 117 
+k_e2 = 339  %%87
+D = 78000  %%88                                                                   
+BW = 78    %%89
+INS_offset = 7 
+be = 3          
+bradykinin = 1  
+pf = 34
+bf_b = 3 
+INS_b = 0                                                               
+p_bf = 1                                                     
+                                            
+V_m = 2.5                                                                       
+V_mx = 0.047                                                                  
+K_m = 225.59                                                               
+k2 = 0.0331 
+k1 = 0.0331           
+part_m = 0.5 
+part_f = 0.44 
+part_l = 1.67 
+U_ii = 0.8447
 
 ********** MODEL VARIABLES
+
+I = I_p/V_I                                                                      
+G = G_p/V_G  
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
-EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_lmax = (1-part)*(V_l0+V_mX*INS)                                                
-V_fmax = part*(V_f0+V_fX*INS)                                                    
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                      
 E = 0                                                                            
-S = gamma*I_po                                                                                                                       
-I = I_p/V_I                                                                      
-G = G_p/V_G                                                                      
+S = gamma*I_po                                                                     
 HE = (-m_5*S)+m_6                                                                
 m_3 = HE*m_1/(1-HE)                                                              
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
 k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)   
-bf_f = (3+0.1*kbf*(INS_f_bf+5)*1e-7)*BRF 
-bfe_f = (bf_f-bf_b)*(INS_f_bf-INS_b)*scale_bf                                                                                                                                                                                 
+bf_f = bradykinin*(be+kbf*(INS_f+INS_offset))
+bfe_f = bf_f %(bf_f-bf_b)*(INS_f-INS_b)*p_bf
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
+                                  
+V_mmax = V_m+V_mx*INS
 
-V_mmax_new = scale4*(V_m0+V_mx_new*INS_new)
-
-U_idm = scale6*V_mmax_new*G_t/(K_m0+G_t)
-U_idf = scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3) 
-U_liver = scale10*V_lmax*G_t/(K_l0+G_t)
+U_idm = part_m*V_mmax*G_t/(K_m+G_t)
+U_idf = part_f*scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3) 
+U_idl = part_l*U_idm 
 
-U_id = U_idf+U_idm+U_liver 
-U = U_ii+U_id
+U_id = U_idf + U_idm + U_idl
 
-U_fm = U_idm+U_idf
+********** MODEL REACTIONS                                                                                                                                
 
-********** MODEL REACTIONS                                                                                                          
-
-v1a      = IR*k1a*(INS+5)*1e-3
+v1a      = IR*k1a*(INS_f+5)*1e-3
 v1basal  = k1basal*IR
 v1c      = IRins*k1c
 v1d      = IRp*k1d
@@ -282,6 +258,9 @@ v9b1     = S6Kp*k9b1
 v9f2     = S6*k9f2*S6Kp
 v9b2     = S6p*k9b2
 
+V_1 = k1*(I-I_b)
+V_2 = k2*INS 
+
 ********** MODEL FUNCTIONS
 
 ********** MODEL EVENTS
diff --git a/scripts/Models/M2b.txt b/scripts/Models/M2b.txt
index 1cf9b1df20cc8444d13df1ccf8c8148538d923e3..a420b7eff5268c33926c90f0bc9cc9a800a7bf5c 100644
--- a/scripts/Models/M2b.txt
+++ b/scripts/Models/M2b.txt
@@ -2,10 +2,9 @@
 M2b
 
 ********** MODEL NOTES
-
+M2b - addition of intracellular glucose phosphorylation to adipocyte module
 
 ********** MODEL STATES
-
 d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
 d/dt(IRp)               = v1basal+v1c-v1d-v1g   %%2
 d/dt(IRins)             = v1a-v1c               %%3
@@ -14,44 +13,42 @@ d/dt(IRi)               = v1e-v1r               %%5
 d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
 d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
 d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
-d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
 d/dt(X)                 = v3b-v3a               %%10
-d/dt(Xp)                = v3a-v3b               %%12
-d/dt(PKB)               = -v4a+v4b+v4h          %%13   
-d/dt(PKB308p)           = v4a-v4b-v4c           %%14
-d/dt(PKB473p)           = -v4e+v4f-v4h          %%15
-d/dt(PKB308p473p)       = v4c+v4e-v4f           %%16
-d/dt(mTORC1)            = v5b-v5a               %%17         
-d/dt(mTORC1a)           = v5a-v5b               %%18
-d/dt(mTORC2)            = -v5c+v5d              %%19
-d/dt(mTORC2a)           = v5c-v5d               %%20
-d/dt(AS160)             = v6b1-v6f1             %%21
-d/dt(AS160p)            = v6f1-v6b1             %%22
-d/dt(GLUT4m)            = (v7f-v7b)             %%23
-d/dt(GLUT4)             = -v7f+v7b              %%24
-d/dt(S6K)               = v9b1-v9f1             %%26
-d/dt(S6Kp)              = v9f1-v9b1             %%27
-d/dt(S6)                = v9b2-v9f2             %%28
-d/dt(S6p)               = v9f2-v9b2             %%29                
+d/dt(Xp)                = v3a-v3b               %%11
+d/dt(PKB)               = -v4a+v4b+v4h          %%12   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%13
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%14
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%15
+d/dt(mTORC1)            = v5b-v5a               %%16         
+d/dt(mTORC1a)           = v5a-v5b               %%17
+d/dt(mTORC2)            = -v5c+v5d              %%18
+d/dt(mTORC2a)           = v5c-v5d               %%19
+d/dt(AS160)             = v6b1-v6f1             %%20
+d/dt(AS160p)            = v6f1-v6b1             %%21
+d/dt(GLUT4m)            = v7f-v7b               %%22
+d/dt(GLUT4)             = -v7f+v7b              %%23
+d/dt(S6K)               = v9b1-v9f1             %%24
+d/dt(S6Kp)              = v9f1-v9b1             %%25
+d/dt(S6)                = v9b2-v9f2             %%26
+d/dt(S6p)               = v9f2-v9b2             %%27             
+
+d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
+d/dt(G6P) = V_G6P-V_met           
                                                       
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t-0.6*(U_liver) 
-d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                           
-d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S                                         
+d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
+d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                      
+d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S  
+d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                            
 d/dt(Q_sto1) = -k_gri*Q_sto1                                                     
 d/dt(Q_sto2) = (-k_empt*Q_sto2)+k_gri*Q_sto1                                     
 d/dt(Q_gut) = (-k_abs*Q_gut)+k_empt*Q_sto2                                       
 d/dt(I_1) = -k_i*(I_1-I)                                                         
-d/dt(I_d) = -k_i*(I_d-I_1)                                                       
-d/dt(INS_p) = (-p_2U*INS_p)+p_2U*(I-I_b)                                             
+d/dt(I_d) = -k_i*(I_d-I_1)                                                           
 d/dt(I_po) = (-gamma*I_po)+S_po                                                  
-d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                
-d/dt(INS) = V_1-V_2
-d/dt(INS_f) = (-C1_f*INS_f)+C2_f*(I-I_b)
-d/dt(INS_f_bf) = (-C1_bf*INS_f_bf)+C2_bf*(I-I_b)
-d/dt(time)     = 1
-d/dt(Glu_in) = scale1*(V_in-V_out)-V_G6P 
-d/dt(G6P) = V_G6P-V_met 
+d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                        
+d/dt(INS_f) = (-p_2U*INS_f)+p_2U*(I-I_b)  
+d/dt(INS) = V_1-V_2    
 
 IR(0)               = 100
 IRp(0)              = 0
@@ -79,7 +76,11 @@ GLUT4(0)            = 100
 S6K(0)              = 100
 S6Kp(0)             = 0
 S6(0)               = 100
-S6p(0)              = 0                               
+S6p(0)              = 0    
+
+Glu_in(0) = 0
+G6P(0) = 0
+                           
 G_p(0) = 178                                                                   
 G_t(0) = 130                                                                     
 I_l(0) = 4.5                                                                 
@@ -89,14 +90,12 @@ Q_sto2(0) = 0
 Q_gut(0) = 0                                                                     
 I_1(0) = 25                                                                      
 I_d(0) = 25                                                                      
-INS_p(0) = 0                                                                     
 I_po(0) = 3.6                                                  
 Y(0) = 0
+INS_f(0) = 0
 INS(0) = 0
-Glu_in(0) = 0
-G6P(0) = 0
 
-********** MODEL PARAMETERS
+********** MODEL PARAMETERS 
 
 diabetes = 1  %%1                                                                
 k1a = 0.633141  %%2                                                              
@@ -142,9 +141,9 @@ k9f2 = 3.3289  %%41
 k9b2 = 30.9967  %%42                                                             
 km9 = 5872.68  %%43                                                              
 n9 = 0.985466  %%44                                                              
-kbf = 1e+06  %%45                                                                
-scaleModel = 2.1e-06  %%46  
-                                                     
+kbf = 0.01 %%45      1e+06                                                           
+nC = 2.1e-06  %%46  
+
 V_G = 1.88  %%47                                                                 
 k_1 = 0.065  %%48                                                                
 k_2 = 0.079  %%49                                                                
@@ -165,19 +164,11 @@ k_gri = 0.0558  %%63
 f = 0.9  %%64                                                                    
 b = 0.82  %%65                                                                   
 d = 0.01  %%66                                                                   
-BW = 78  %%67                                                                    
 k_p1 = 2.7  %%68                                                                 
 k_p2 = 0.0021  %%69                                                              
 k_p3 = 0.009  %%70                                                               
 k_p4 = 0.0618  %%71                                                              
 k_i = 0.0079  %%72                                                               
-U_ii = 1.23  %%73                                                                
-V_l0 = 2.5     %% 74                                                                  
-V_mX = 0.047  %%75                                                               
-K_l0 = 225.59  %%76                                                              
-V_f0 = 2.5  %%77                                                                 
-V_fX = 0.047  %%78                                                               
-K_f0 = 225.59  %%79                                                              
 p_2U = 0.0331  %%80                                                              
 part = 0.2  %%81                                                                 
 K = 2.3  %%82                                                                    
@@ -185,85 +176,74 @@ alpha = 0.05  %%83
 beta = 0.11  %%84                                                                
 gamma = 0.5  %%85                                                                
 k_e1 = 0.0005  %%86                                                              
-k_e2 = 339  %%87                                                                 
-D = 78000  %%88    
-
-scale10 = 0.5 %% 89
-Ip_b = 7 %% 90
-be = 3 %% 91
-
-k2 = 2.15 %% 90
-k1 = 1.3 %% 91
-C1_f = 1.25 %% 92
-C2_f = 0.6 %% 93
-V_m0 = 2.5 %% 94
-V_mx_new = 0.047 %% 96
-K_m0 = 225.59 %% 98
-BRF = 1 %% 100
-bf_b = 3 %% 101
-INS_b = 0.8549 %% 102
-scale_bf = 0.05e2 %% 103
-C1_bf = 0.075 %% 104
-C2_bf = 0.06 %% 105
-scale1 = 1 %% 106
-scale2 = 0.13 %% 107
-scale3 = 2.4 %% 108
-scale4 = 0.4 %% 109
-scale6 = 1.5 %% 110
-scale7 = 0.7 %% 111
-scale8 = 1.8 %% 112
-scale9 = 14 %% 113
-k_gluin = 1300 %% 114
-k_G6P = 15 %% 115
-V_G6Pmax = 90 %% 116
-pl = 0.5 %% 117  
+k_e2 = 339  %%87
+D = 78000  %%88                                                                   
+BW = 78    %%89
+INS_offset = 7 
+be = 3          
+pf = 34
+bf_b = 3 
+INS_b = 0                                                               
+p_bf = 1                                                     
+                                            
+U_ii = 1.47  
+k2 = 0.0331 
+k1 = 0.0331
+V_m = 2.5 
+V_mx = 0.047 
+K_m = 225.59 
+V_l = 2.5 
+V_lx = 0.047 
+K_l = 225.59
+p1 = 1 
+p2 = 2.4 
+p3 = 0.7
+p4 = 1.8  
+k_gluin = 1300
+k_G6P = 15 
+V_G6Pmax = 90  
 
 ********** MODEL VARIABLES
+
+I = I_p/V_I                                                                      
+G = G_p/V_G  
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
-EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_lmax = (1-part)*(V_l0+V_mX*G_t)                                                
-V_fmax = part*(V_f0+V_fX*INS_p)                                                    
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                      
 E = 0                                                                            
-S = gamma*I_po 
-I = I_p/V_I                                                                      
-G = G_p/V_G                                                                      
+S = gamma*I_po                                                                     
 HE = (-m_5*S)+m_6                                                                
 m_3 = HE*m_1/(1-HE)                                                              
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
 k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)   
-bf_f = (3+0.1*kbf*(INS_f_bf+5)*1e-7)*BRF 
-bfe_f = (bf_f-bf_b)*(INS_f_bf-INS_b)*scale_bf
+bf_f = (be+kbf*(INS_f+INS_offset))
+bfe_f = bf_f
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
+                                  
+V_mmax = V_m+V_mx*INS
+V_lmax = V_l+V_lx*INS
 
-V_mmax = scale4*(V_m0+V_mx_new*INS)
-
-U_idf_cell = scaleModel*(k8*GLUT4m/34 + glut1/34+ bfe_f)
-V_in = scale8*(G_t-Glu_in)*U_idf_cell
-V_out = scale7*Glu_in
-
-U_idm = scale6*V_mmax*G_t/(K_m0+G_t)
-U_idf = scale9*(V_in-V_out)
-U_liver = scale10*V_lmax*G_t/(K_l0+G_t)
+INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
+V_in = p4*G_t*INS_fe
+V_out = p3*Glu_in
 
-U_id = U_idf+U_idm 
-U = U_ii+U_id+U_liver
+U_idm = V_mmax*G_t/(K_m+G_t)
+U_idl = V_lmax*G_t/(K_l+G_t)
+U_idf = V_in-V_out 
 
-yf=U_idf
-ym=U_idm
+U_id = U_idf + U_idm + U_idl
 
 ********** MODEL REACTIONS                                                                                                                                
 
-v1a      = IR*k1a*(INS_p+5)*1e-3
+v1a      = IR*k1a*(INS_f+5)*1e-3
 v1basal  = k1basal*IR
 v1c      = IRins*k1c
 v1d      = IRp*k1d
 v1e      = IRip*k1f*Xp
 v1g      = IRp*k1g
 v1r      = IRi*k1r
-v2a      = IRS1*k2a*
-IRip
+v2a      = IRS1*k2a*IRip
 v2b      = IRS1p*k2b
 v2c      = IRS1p*k2c*mTORC1a*diabetes
 v2d      = IRS1p307*k2d
@@ -295,13 +275,13 @@ v9f2     = S6*k9f2*S6Kp
 v9b2     = S6p*k9b2
 
 V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
-V_met = scale3*G6P
+V_met = p2*G6P
 
 V_1 = k1*(I-I_b)
-V_2 = k2*INS
+V_2 = k2*INS 
 
 ********** MODEL FUNCTIONS
 
 ********** MODEL EVENTS
 
-********** MODEL MATLAB FUNCTIONS
+********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/suport/IQMtools/IQMlite/classeshandling/modelhandling/auxiliary/SBMLexport/M1.txt b/scripts/Models/M3.txt
similarity index 75%
rename from scripts/suport/IQMtools/IQMlite/classeshandling/modelhandling/auxiliary/SBMLexport/M1.txt
rename to scripts/Models/M3.txt
index eff9181393930c4444c7c6a339ce0a6bb0c00ac4..dc131d11fdafd72f28b996cd97e1d25993bb57e8 100644
--- a/scripts/suport/IQMtools/IQMlite/classeshandling/modelhandling/auxiliary/SBMLexport/M1.txt
+++ b/scripts/Models/M3.txt
@@ -1,8 +1,8 @@
 ********** MODEL NAME
-M1
+M3
 
 ********** MODEL NOTES
-
+M3 - blood flow
 
 ********** MODEL STATES
 d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
@@ -13,39 +13,42 @@ d/dt(IRi)               = v1e-v1r               %%5
 d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
 d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
 d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
-d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
 d/dt(X)                 = v3b-v3a               %%10
-d/dt(Xp)                = v3a-v3b               %%12
-d/dt(PKB)               = -v4a+v4b+v4h          %%13   
-d/dt(PKB308p)           = v4a-v4b-v4c           %%14
-d/dt(PKB473p)           = -v4e+v4f-v4h          %%15
-d/dt(PKB308p473p)       = v4c+v4e-v4f           %%16
-d/dt(mTORC1)            = v5b-v5a               %%17         
-d/dt(mTORC1a)           = v5a-v5b               %%18
-d/dt(mTORC2)            = -v5c+v5d              %%19
-d/dt(mTORC2a)           = v5c-v5d               %%20
-d/dt(AS160)             = v6b1-v6f1             %%21
-d/dt(AS160p)            = v6f1-v6b1             %%22
-d/dt(GLUT4m)            = (v7f-v7b)             %%23
-d/dt(GLUT4)             = -v7f+v7b              %%24
-d/dt(S6K)               = v9b1-v9f1             %%26
-d/dt(S6Kp)              = v9f1-v9b1             %%27
-d/dt(S6)                = v9b2-v9f2             %%28
-d/dt(S6p)               = v9f2-v9b2             %%29                                               
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                        
-d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S                                         
-d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                         
+d/dt(Xp)                = v3a-v3b               %%11
+d/dt(PKB)               = -v4a+v4b+v4h          %%12   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%13
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%14
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%15
+d/dt(mTORC1)            = v5b-v5a               %%16         
+d/dt(mTORC1a)           = v5a-v5b               %%17
+d/dt(mTORC2)            = -v5c+v5d              %%18
+d/dt(mTORC2a)           = v5c-v5d               %%19
+d/dt(AS160)             = v6b1-v6f1             %%20
+d/dt(AS160p)            = v6f1-v6b1             %%21
+d/dt(GLUT4m)            = v7f-v7b               %%22
+d/dt(GLUT4)             = -v7f+v7b              %%23
+d/dt(S6K)               = v9b1-v9f1             %%24
+d/dt(S6Kp)              = v9f1-v9b1             %%25
+d/dt(S6)                = v9b2-v9f2             %%26
+d/dt(S6p)               = v9f2-v9b2             %%27             
+
+d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
+d/dt(G6P) = V_G6P-V_met           
+                                                      
+d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
+d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                      
+d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S  
+d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                            
 d/dt(Q_sto1) = -k_gri*Q_sto1                                                     
 d/dt(Q_sto2) = (-k_empt*Q_sto2)+k_gri*Q_sto1                                     
 d/dt(Q_gut) = (-k_abs*Q_gut)+k_empt*Q_sto2                                       
 d/dt(I_1) = -k_i*(I_1-I)                                                         
-d/dt(I_d) = -k_i*(I_d-I_1)                                                       
-d/dt(INS) = (-p_2U*INS)+p_2U*(I-I_b)                                             
+d/dt(I_d) = -k_i*(I_d-I_1)                                                           
 d/dt(I_po) = (-gamma*I_po)+S_po                                                  
-d/dt(Y) = -alpha*(Y-beta*(G-G_b))  
-d/dt(E) = 0                                              
-
+d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                        
+d/dt(INS_f) = (-p_2U*INS_f)+p_2U*(I-I_b)  
+d/dt(INS) = V_1-V_2    
 
 IR(0)               = 100
 IRp(0)              = 0
@@ -73,7 +76,11 @@ GLUT4(0)            = 100
 S6K(0)              = 100
 S6Kp(0)             = 0
 S6(0)               = 100
-S6p(0)              = 0                               
+S6p(0)              = 0    
+
+Glu_in(0) = 0
+G6P(0) = 0
+                           
 G_p(0) = 178                                                                   
 G_t(0) = 130                                                                     
 I_l(0) = 4.5                                                                 
@@ -83,12 +90,12 @@ Q_sto2(0) = 0
 Q_gut(0) = 0                                                                     
 I_1(0) = 25                                                                      
 I_d(0) = 25                                                                      
-INS(0) = 0                                                                     
 I_po(0) = 3.6                                                  
 Y(0) = 0
-E(0) = 0
+%INS_f(0) = 0
+INS(0) = 0
 
-********** MODEL PARAMETERS
+********** MODEL PARAMETERS 
 
 diabetes = 1  %%1                                                                
 k1a = 0.633141  %%2                                                              
@@ -134,9 +141,16 @@ k9f2 = 3.3289  %%41
 k9b2 = 30.9967  %%42                                                             
 km9 = 5872.68  %%43                                                              
 n9 = 0.985466  %%44                                                              
-kbf = 1e+06  %%45                                                                
-scaleModel = 2.1e-06  %%46  
-                                                     
+kbf = 0.01 %%45      1e+06                                                           
+nC = 2.1e-06  %%46  
+
+INS_offset = 7 
+be = 3          
+bradykinin = 1  
+pf = 34
+bf_b = 3 
+INS_b = 0
+
 V_G = 1.88  %%47                                                                 
 k_1 = 0.065  %%48                                                                
 k_2 = 0.079  %%49                                                                
@@ -157,19 +171,11 @@ k_gri = 0.0558  %%63
 f = 0.9  %%64                                                                    
 b = 0.82  %%65                                                                   
 d = 0.01  %%66                                                                   
-BW = 78  %%67                                                                    
 k_p1 = 2.7  %%68                                                                 
 k_p2 = 0.0021  %%69                                                              
 k_p3 = 0.009  %%70                                                               
 k_p4 = 0.0618  %%71                                                              
 k_i = 0.0079  %%72                                                               
-
-V_m0 = 2.5    %%74                                                                   
-V_mX = 0.047  %%75                                                                
-K_m0 = 225.59  %%76                                             
-V_f0 = 2.5  %%77                                                                 
-V_fX = 0.047  %%78                                                               
-K_f0 = 225.59  %%79                                                              
 p_2U = 0.0331  %%80                                                              
 part = 0.2  %%81                                                                 
 K = 2.3  %%82                                                                    
@@ -177,39 +183,62 @@ alpha = 0.05  %%83
 beta = 0.11  %%84                                                                
 gamma = 0.5  %%85                                                                
 k_e1 = 0.0005  %%86                                                              
-k_e2 = 339  %%87                                                                 
-D = 78000  %%88    
-
-K_l0 = 225.59
-V_l0 = 2.5                                                                       
-V_lX = 0.047  
+k_e2 = 339  %%87
+D = 78000  %%88                                                                   
+BW = 78                                                                   
+p_bf = 1                                                     
+                                            
+U_ii = 1.47  
+k2 = 0.0331 
+k1 = 0.0331
+V_m = 2.5 
+V_mx = 0.047 
+K_m = 225.59 
+V_l = 2.5 
+V_lx = 0.047 
+K_l = 225.59
+p1 = 1 
+p2 = 2.4 
+p3 = 0.7
+p4 = 1.8  
+k_gluin = 1300
+k_G6P = 15 
+V_G6Pmax = 90  
 
 ********** MODEL VARIABLES
+
+I = I_p/V_I                                                                      
+G = G_p/V_G  
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
-EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_mmax = (1-part)*(V_m0+V_mX*INS)                                                
-V_fmax = part*(V_f0+V_fX*INS) 
-V_lmax = (1-part)*(V_l0+V_lX*G_t)                                                   
-Ee = k_e1*(G_p - k_e2)                                                                            
-S = gamma*I_po                                                                                                                       
-I = I_p/V_I                                                                      
-G = G_p/V_G                                                                      
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                      
+E = 0                                                                            
+S = gamma*I_po                                                                     
 HE = (-m_5*S)+m_6                                                                
 m_3 = HE*m_1/(1-HE)                                                              
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
-k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)     
-U_idm = 0.5*V_mmax*G_t/(K_m0+G_t)
-U_idf = 0.5*scaleModel*(k8*GLUT4m*G_t/34 + glut1*G_t/34 + kbf*(INS+5)*1e-3)                                                                                                              
-U = U_ii+U_id                                                                    
+k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)   
+bf_f = bradykinin*(be+kbf*(INS_f+INS_offset))
+bfe_f = (bf_f-bf_b)*(INS_f-INS_b)*p_bf
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
-U_liver = 45*U_idm/27
-U_ii = 23*U_liver/45
-U_id = U_idm+U_idf+U_liver
+                                  
+V_mmax = V_m+V_mx*INS
+V_lmax = V_l+V_lx*INS
+
+INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
+V_in = p4*G_t*INS_fe
+V_out = p3*Glu_in
+
+U_idm = V_mmax*G_t/(K_m+G_t)
+U_idl = V_lmax*G_t/(K_l+G_t)
+U_idf = V_in-V_out 
 
-********** MODEL REACTIONS                                                                                                          
-v1a      = IR*k1a*(INS+5)*1e-3
+U_id = U_idf + U_idm + U_idl
+
+********** MODEL REACTIONS                                                                                                                                
+
+v1a      = IR*k1a*(INS_f+5)*1e-3
 v1basal  = k1basal*IR
 v1c      = IRins*k1c
 v1d      = IRp*k1d
@@ -247,11 +276,14 @@ v9b1     = S6Kp*k9b1
 v9f2     = S6*k9f2*S6Kp
 v9b2     = S6p*k9b2
 
-glucoseuptake = k8*GLUT4m*(G_t*5/170) + glut1*(G_t*5/170) + kbf*((INS+5)*1e-3)
+V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
+V_met = p2*G6P
+
+V_1 = k1*(I-I_b)
+V_2 = k2*INS 
 
 ********** MODEL FUNCTIONS
 
 ********** MODEL EVENTS
-event1 = gt(G_p,339),E,Ee
-event2 = lt(G_p,339),E,0
+
 ********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/Models/M3m.txt b/scripts/Models/M3m.txt
deleted file mode 100644
index a644fec5a77e880a7d59a846b1218682ed514132..0000000000000000000000000000000000000000
--- a/scripts/Models/M3m.txt
+++ /dev/null
@@ -1,128 +0,0 @@
-********** MODEL NAME
-M3m
-
-********** MODEL NOTES
-sub-model version of M3, with EGP, Ra, GLUT4m, and INS as input curves
-
-********** MODEL STATES
-          
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t    
-d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
-d/dt(G6P) = V_G6P-V_met   
-
-G_p(0) = 178                                                                   
-G_t(0) = 130  
-Glu_in(0) = 0
-G6P(0) = 0
-
-********** MODEL PARAMETERS
-            
-bradykinin = 1   
-kbf = 0.01
-INS_offset = 7 
-be = 3                                                            
-nC = 2.1e-06 
-k8 = 724.242                                        
-glut1 = 7042.19
-k_1 = 0.065
-k_2 = 0.079 
-pf = 34
-INS_b = 0
-p_bf = 1                                      
-bf_b = 3 
-
-EGP1 = 1
-EGP2 = 1
-EGP3 = 1
-EGP4 = 1
-EGP5 = 1
-EGP6 = 1
-EGP7 = 1
-EGP8 = 1
-EGP9 = 1
-
-Ra1 = 1
-Ra2 = 1
-Ra3 = 1
-Ra4 = 1
-Ra5 = 1
-Ra6 = 1
-Ra7 = 1
-Ra8 = 1
-Ra9 = 1
-
-INS1 = 1
-INS2 = 1
-INS3 = 1
-INS4 = 1
-INS5 = 1
-INS6 = 1
-INS7 = 1
-INS8 = 1
-INS9 = 1
-
-GLm1 = 1
-GLm2 = 1 
-GLm3 = 1 
-GLm4 = 1 
-GLm5 = 1 
-GLm6 = 1 
-GLm7 = 1 
-GLm8 = 1 
-GLm9 = 1
- 
-V_m = 2.5 
-V_mx = 0.047 
-K_m = 225.59 
-part_l = 0.5
-K_l = 225.59
-p1 = 1 
-p2 = 2.4 
-p3 = 0.7
-p4 = 1.8 
-part_f = 14  
-k_gluin = 1300
-k_G6P = 15 
-V_G6Pmax = 90 
-U_ii = 1.23
-
-********** MODEL VARIABLES
-
-INS = interpcsIQM([0,30,60,90,120,180,240,300,360],[INS1,INS2,INS3,INS4,INS5,INS6,INS7,INS8,INS9],time)
-EGP = interpcsIQM([0,30,60,90,120,180,240,300,360],[EGP1,EGP2,EGP3,EGP4,EGP5,EGP6,EGP7,EGP8,EGP9],time)
-Ra = interpcsIQM([0,30,60,90,120,180,240,300,360],[Ra1,Ra2,Ra3,Ra4,Ra5,Ra6,Ra7,Ra8,Ra9],time)
-GLUT4m = interpcsIQM([0,30,60,90,120,180,240,300,360],[GLm1,GLm2,GLm3,GLm4,GLm5,GLm6,GLm7,GLm8,GLm9],time)
-
-E = 0
-
-bf_f = bradykinin*(be+kbf*(INS+INS_offset))
-bfe_f = (bf_f-bf_b)*(INS-INS_b)*p_bf  
-
-part_m = 1-part_l                                     
-V_mmax = part_m*(V_m+V_mx*INS)
-V_lmax = part_l*V_mmax %(V_l+V_lx*INS) 
-
-INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
-V_in = p4*G_t*INS_fe
-V_out = p3*Glu_in
-
-U_idm = V_mmax*G_t/(K_m+G_t)
-U_idl = V_lmax*G_t/(K_l+G_t)
-U_idf = part_f*(V_in-V_out) 
-
-U_id = U_idf + U_idm + U_idl
-
-V_G = 1.88
-G = G_p/V_G
-
-********** MODEL REACTIONS                                                                                                                                
-
-V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
-V_met = p2*G6P
-
-********** MODEL FUNCTIONS
-
-********** MODEL EVENTS
-
-********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/Models/M4.txt b/scripts/Models/M4.txt
index ea0df87f690abdc3a8c4295529095963457740cd..25b2b22bfc69c124c8f160f2bba72cb1f32ada25 100644
--- a/scripts/Models/M4.txt
+++ b/scripts/Models/M4.txt
@@ -2,10 +2,9 @@
 M4
 
 ********** MODEL NOTES
-
+M4 - final model
 
 ********** MODEL STATES
-
 d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
 d/dt(IRp)               = v1basal+v1c-v1d-v1g   %%2
 d/dt(IRins)             = v1a-v1c               %%3
@@ -14,7 +13,7 @@ d/dt(IRi)               = v1e-v1r               %%5
 d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
 d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
 d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
-d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
 d/dt(X)                 = v3b-v3a               %%10
 d/dt(Xp)                = v3a-v3b               %%11
 d/dt(PKB)               = -v4a+v4b+v4h          %%12   
@@ -27,29 +26,29 @@ d/dt(mTORC2)            = -v5c+v5d              %%18
 d/dt(mTORC2a)           = v5c-v5d               %%19
 d/dt(AS160)             = v6b1-v6f1             %%20
 d/dt(AS160p)            = v6f1-v6b1             %%21
-d/dt(GLUT4m)            = (v7f-v7b)             %%22
+d/dt(GLUT4m)            = v7f-v7b               %%22
 d/dt(GLUT4)             = -v7f+v7b              %%23
 d/dt(S6K)               = v9b1-v9f1             %%24
 d/dt(S6Kp)              = v9f1-v9b1             %%25
 d/dt(S6)                = v9b2-v9f2             %%26
-d/dt(S6p)               = v9f2-v9b2             %%27                
-                 
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t                                        
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t-U_idl                                        
-d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S                                         
-d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                         
+d/dt(S6p)               = v9f2-v9b2             %%27             
+
+d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
+d/dt(G6P) = V_G6P-V_met           
+                                                      
+d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
+d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                      
+d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S  
+d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                            
 d/dt(Q_sto1) = -k_gri*Q_sto1                                                     
 d/dt(Q_sto2) = (-k_empt*Q_sto2)+k_gri*Q_sto1                                     
 d/dt(Q_gut) = (-k_abs*Q_gut)+k_empt*Q_sto2                                       
 d/dt(I_1) = -k_i*(I_1-I)                                                         
-d/dt(I_d) = -k_i*(I_d-I_1)                                                       
-d/dt(INS_p) = (-p_2U*INS_p)+p_2U*(I-I_b)                                             
+d/dt(I_d) = -k_i*(I_d-I_1)                                                           
 d/dt(I_po) = (-gamma*I_po)+S_po                                                  
-d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                
-d/dt(INS) = V_1-V_2
-d/dt(INS_f_bf) = (-C1_bf*INS_f_bf)+C2_bf*(I-I_b)
-d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
-d/dt(G6P) = V_G6P-V_met 
+d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                        
+d/dt(INS_f) = (-p_2U*INS_f)+p_2U*(I-I_b)  
+d/dt(INS) = V_1-V_2    
 
 IR(0)               = 100
 IRp(0)              = 0
@@ -77,7 +76,11 @@ GLUT4(0)            = 100
 S6K(0)              = 100
 S6Kp(0)             = 0
 S6(0)               = 100
-S6p(0)              = 0                               
+S6p(0)              = 0    
+
+Glu_in(0) = 0
+G6P(0) = 0
+                           
 G_p(0) = 178                                                                   
 G_t(0) = 130                                                                     
 I_l(0) = 4.5                                                                 
@@ -87,15 +90,12 @@ Q_sto2(0) = 0
 Q_gut(0) = 0                                                                     
 I_1(0) = 25                                                                      
 I_d(0) = 25                                                                      
-INS_p(0) = 0                                                                     
 I_po(0) = 3.6                                                  
 Y(0) = 0
+INS_f(0) = 0
 INS(0) = 0
-INS_f_bf(0) = 0
-Glu_in(0) = 0
-G6P(0) = 0
 
-********** MODEL PARAMETERS
+********** MODEL PARAMETERS 
 
 diabetes = 1  %%1                                                                
 k1a = 0.633141  %%2                                                              
@@ -141,9 +141,9 @@ k9f2 = 3.3289  %%41
 k9b2 = 30.9967  %%42                                                             
 km9 = 5872.68  %%43                                                              
 n9 = 0.985466  %%44                                                              
-kbf = 0.01  %%45                                                                
+kbf = 0.01 %%45      1e+06                                                           
 nC = 2.1e-06  %%46  
-                                                     
+
 V_G = 1.88  %%47                                                                 
 k_1 = 0.065  %%48                                                                
 k_2 = 0.079  %%49                                                                
@@ -164,95 +164,80 @@ k_gri = 0.0558  %%63
 f = 0.9  %%64                                                                    
 b = 0.82  %%65                                                                   
 d = 0.01  %%66                                                                   
-BW = 78  %%67                                                                    
-k_p1 = 2.7  %%68                                                                 
-k_p2 = 0.0021  %%69                                                              
-k_p3 = 0.009  %%70                                                               
-k_p4 = 0.0618  %%71                                                              
-k_i = 0.0079  %%72                                                               
-U_ii = 1.23  %%73                                                                
-V_l = 2.5     %% 74                                                                  
-V_lX = 0.047  %%75                                                               
-K_l = 225.59  %%76                                                              
-V_f = 2.5  %%77                                                                 
-V_fX = 0.047  %%78                                                               
-K_f0 = 225.59  %%79                                                              
-p_2U = 0.0331  %%80                                                              
-part = 0.2  %%81                                                                 
-K = 2.3  %%82                                                                    
-alpha = 0.05  %%83                                                               
-beta = 0.11  %%84                                                                
-gamma = 0.5  %%85                                                                
-k_e1 = 0.0005  %%86                                                              
-k_e2 = 339  %%87                                                                 
-D = 78000  %%88    
+k_p1 = 2.7  %%67                                                                 
+k_p2 = 0.0021  %%68                                                              
+k_p3 = 0.009  %%69                                                               
+k_p4 = 0.0618  %%70                                                              
+k_i = 0.0079  %%71                                                      
+p_2U = 0.0331  %%72                                                               
+K = 2.3  %%73                                                                    
+alpha = 0.05  %%74                                                               
+beta = 0.11  %%75                                                                
+gamma = 0.5  %%76                                                                
+k_e1 = 0.0005  %%77                                                              
+k_e2 = 339  %%78
+D = 78000  %%79                                                                   
+BW = 78    %%80
 
-p10 = 0.5 %% 89
-INS_offset = 7 %% 90
-be = 3 %% 91
-
-k2 = 2.15 %% 92
-k1 = 1.3 %% 93
-C1_f = 1.25 %% 94
-C2_f = 0.6 %% 95
-V_m = 2.5 %% 96
-V_mx_new = 0.047 %% 97
-K_m = 225.59 %% 98
-bradykinin = 1 %% 99
-bf_b = 3 %% 100
-INS_b = 0.8549 %% 101
-p_bf = 0.05e2 %% 102
-C1_bf = 0.075 %% 103
-C2_bf = 0.06 %% 104
-p1 = 1 %% 105
-p2 = 0.13 %% 106
-p3 = 2.4 %% 107
-p4 = 0.4 %% 108
-p6 = 1.5 %% 109
-p7 = 0.7 %% 110
-p8 = 1.8 %% 111
-p9 = 14 %% 112
-k_gluin = 1300 %% 113
-k_G6P = 15 %% 114
-V_G6Pmax = 90 %% 115
-pl = 0.5 %% 116
-pf = 34 %% 117
+INS_offset = 7 
+be = 3          
+bradykinin = 1  
+pf = 34
+bf_b = 3 
+INS_b = 0                                                               
+p_bf = 1                                                     
+                                            
+U_ii = 1.47  
+k2 = 0.0331 
+k1 = 0.0331
+V_m = 2.5 
+V_mx = 0.047 
+K_m = 225.59 
+V_l = 2.5 
+V_lx = 0.047 
+K_l = 225.59
+p1 = 1 
+p2 = 2.4 
+p3 = 0.7
+p4 = 1.8  
+k_gluin = 1300
+k_G6P = 15 
+V_G6Pmax = 90  
 
 ********** MODEL VARIABLES
+
+I = I_p/V_I                                                                      
+G = G_p/V_G  
 aa = 5/2/(1-b)/D                                                                 
 cc = 5/2/d/D                                                                     
-EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                           
-V_lmax = (1-part)*(V_l+V_lX*INS_p)                                                
-V_fmax = part*(V_f+V_fX*INS_p)                                                    
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                      
 E = 0                                                                            
-S = gamma*I_po 
-I = I_p/V_I                                                                      
-G = G_p/V_G                                                                      
+S = gamma*I_po                                                                     
 HE = (-m_5*S)+m_6                                                                
 m_3 = HE*m_1/(1-HE)                                                              
 Q_sto = Q_sto1+Q_sto2                                                            
 Ra = f*k_abs*Q_gut/BW                                                           
 k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)   
-bf_f = bradykinin*(be+kbf*(INS_f_bf+INS_offset))
-bfe_f = (bf_f-bf_b)*(INS_f_bf-INS_b)*p_bf
+bf_f = bradykinin*(be+kbf*(INS_f+INS_offset))
+bfe_f = bf_f %(bf_f-bf_b)*(INS_f-INS_b)*p_bf
 S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
-
-V_mmax = p4*(V_m+V_mx_new*INS)
+                                  
+V_mmax = V_m+V_mx*INS
+V_lmax = V_l+V_lx*INS
 
 INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
-V_in = p8*(G_t-Glu_in)*INS_fe
-V_out = p7*Glu_in
-
-U_idm = p6*V_mmax*G_t/(K_m+G_t)
-U_idf = p9*(V_in-V_out)
-U_idl = pl*V_lmax*G_t/(K_l+G_t)
+V_in = p4*G_t*INS_fe
+V_out = p3*Glu_in
 
-U_id = U_idf+U_idm 
+U_idm = V_mmax*G_t/(K_m+G_t)
+U_idl = V_lmax*G_t/(K_l+G_t)
+U_idf = V_in-V_out 
 
+U_id = U_idf + U_idm + U_idl
 
 ********** MODEL REACTIONS                                                                                                                                
 
-v1a      = IR*k1a*(INS_p+5)*1e-3
+v1a      = IR*k1a*(INS_f+5)*1e-3
 v1basal  = k1basal*IR
 v1c      = IRins*k1c
 v1d      = IRp*k1d
@@ -291,13 +276,13 @@ v9f2     = S6*k9f2*S6Kp
 v9b2     = S6p*k9b2
 
 V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
-V_met = p3*G6P
+V_met = p2*G6P
 
 V_1 = k1*(I-I_b)
-V_2 = k2*INS
+V_2 = k2*INS 
 
 ********** MODEL FUNCTIONS
 
 ********** MODEL EVENTS
 
-********** MODEL MATLAB FUNCTIONS
+********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/Models/M1.xml b/scripts/Models/M4_SBML.xml
similarity index 89%
rename from scripts/Models/M1.xml
rename to scripts/Models/M4_SBML.xml
index b7e1f97f84fb19f7c9537b09ba60fcc762b16fba..1663cce653e9fd5a9df4be656507247b98608a63 100644
--- a/scripts/Models/M1.xml
+++ b/scripts/Models/M4_SBML.xml
@@ -1,24 +1,17 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <sbml xmlns="http://www.sbml.org/sbml/level2" level="2" version="1">
-  <model name="M1_SBML">
-    <notes>
-      <html xmlns="http://www.w3.org/1999/xhtml">
-        <body>
-          <p>adipocyte = 1 {isCompartment:}</p>
-          <p>liver = 1 {isCompartment:}</p>
-          <p>tissue = 1 {isCompartment:}</p>
-          <p>muscle = 1 {isCompartment:}</p>
-          <p>plasma = 1 {isCompartment:}</p>
-          <p>fat = 1 {isCompartment:}</p>
-          <p>stomach = 1 {isCompartment:}</p>
-          <p>intestine = 1 {isCompartment:}</p>
-          <p>betaCell = 1 {isCompartment:}</p>
-          <p>kidney = 1 {isCompartment:}</p>
-        </body>
-      </html>
-    </notes>
+  <model name="M4_SBML">
     <listOfCompartments>
-      <compartment id="rootCompartment" name="rootCompartment" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="adipocyte" name="adipocyte" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="liver" name="liver" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="tissue" name="tissue" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="muscle" name="muscle" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="plasma" name="plasma" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="fat" name="fat" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="stomach" name="stomach" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="intestine" name="intestine" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="betaCell" name="betaCell" spatialDimensions="3" size="1" constant="true"/>
+      <compartment id="kidney" name="kidney" spatialDimensions="3" size="1" constant="true"/>
     </listOfCompartments>
     <listOfSpecies>
       <species id="IR" name="IR" compartment="adipocyte" initialConcentration="100" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false">
@@ -237,27 +230,28 @@
           </html>
         </notes>
       </species>
+      <species id="Glu_in" name="Glu_in" compartment="adipocyte" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
+      <species id="G6P" name="G6P" compartment="adipocyte" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
       <species id="G_p" name="G_p" compartment="plasma" initialAmount="178" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="G_t" name="G_t" compartment="tissue" initialAmount="130" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I_l" name="I_l" compartment="liver" initialAmount="4.5" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I_p" name="I_p" compartment="plasma" initialAmount="1.25" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="Q_sto1" name="Q_sto1" compartment="stomach" initialAmount="78000" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="Q_sto2" name="Q_sto2" compartment="stomach" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="Q_gut" name="Q_gut" compartment="intenstine" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
+      <species id="Q_gut" name="Q_gut" compartment="intestine" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I_1" name="I_1" compartment="liver" initialAmount="25" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I_d" name="I_d" compartment="liver" initialAmount="25" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="INS" name="INS" compartment="tissue" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I_po" name="I_po" compartment="liver" initialAmount="3.6" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="Y" name="Y" compartment="betaCell" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
+      <species id="INS_f" name="INS_f" compartment="tissue" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
+      <species id="INS" name="INS" compartment="tissue" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="I" name="I" compartment="plasma" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
       <species id="G" name="G" compartment="plasma" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
       <species id="Q_sto" name="Q_sto" compartment="stomach" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
       <species id="U_idm" name="U_idm" compartment="muscle" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_idf" name="U_idf" compartment="fat" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="S_po" name="S_po" compartment="liver" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_idl" name="U_idl" compartment="liver" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
+      <species id="U_idl" name="U_idl" compartment="fat" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
+      <species id="U_idf" name="U_idf" compartment="liver" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
       <species id="U_id" name="U_id" compartment="tissue" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="glucoseuptake" name="glucoseuptake" compartment="adipocyte" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
     </listOfSpecies>
     <listOfParameters>
       <parameter id="diabetes" name="diabetes" value="1" constant="true">
@@ -612,7 +606,7 @@
           </html>
         </notes>
       </parameter>
-      <parameter id="kbf" name="kbf" value="1000000" constant="true">
+      <parameter id="kbf" name="kbf" value="0.01" constant="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
@@ -632,7 +626,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 47</body>
+%47</body>
           </html>
         </notes>
       </parameter>
@@ -640,7 +634,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 48</body>
+%48</body>
           </html>
         </notes>
       </parameter>
@@ -648,7 +642,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 49</body>
+%49</body>
           </html>
         </notes>
       </parameter>
@@ -656,7 +650,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 50</body>
+%50</body>
           </html>
         </notes>
       </parameter>
@@ -664,7 +658,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 51</body>
+%51</body>
           </html>
         </notes>
       </parameter>
@@ -672,7 +666,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 52</body>
+%52</body>
           </html>
         </notes>
       </parameter>
@@ -696,7 +690,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 55</body>
+%55</body>
           </html>
         </notes>
       </parameter>
@@ -712,7 +706,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 57</body>
+%57</body>
           </html>
         </notes>
       </parameter>
@@ -720,7 +714,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 58</body>
+%58</body>
           </html>
         </notes>
       </parameter>
@@ -728,7 +722,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 59</body>
+%59</body>
           </html>
         </notes>
       </parameter>
@@ -736,7 +730,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 60</body>
+%60</body>
           </html>
         </notes>
       </parameter>
@@ -744,7 +738,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 61</body>
+%61</body>
           </html>
         </notes>
       </parameter>
@@ -752,7 +746,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 62</body>
+%62</body>
           </html>
         </notes>
       </parameter>
@@ -760,7 +754,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 63</body>
+%63</body>
           </html>
         </notes>
       </parameter>
@@ -768,7 +762,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 64</body>
+%64</body>
           </html>
         </notes>
       </parameter>
@@ -776,7 +770,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 65</body>
+%65</body>
           </html>
         </notes>
       </parameter>
@@ -784,15 +778,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 66</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="BW" name="BW" value="78" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-% 67</body>
+%66</body>
           </html>
         </notes>
       </parameter>
@@ -800,7 +786,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 68</body>
+%67</body>
           </html>
         </notes>
       </parameter>
@@ -808,7 +794,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 69</body>
+%68</body>
           </html>
         </notes>
       </parameter>
@@ -816,7 +802,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 70</body>
+%69</body>
           </html>
         </notes>
       </parameter>
@@ -824,7 +810,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 71</body>
+%70</body>
           </html>
         </notes>
       </parameter>
@@ -832,55 +818,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-% 72</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="V_m0" name="V_m0" value="2.5" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%73</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="V_mX" name="V_mX" value="0.047" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%74</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="K_m0" name="K_m0" value="225.59" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%75</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="V_f0" name="V_f0" value="2.5" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%76</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="V_fX" name="V_fX" value="0.047" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%77</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="K_f0" name="K_f0" value="225.59" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%78</body>
+%71</body>
           </html>
         </notes>
       </parameter>
@@ -888,15 +826,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%79</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="part" name="part" value="0.2" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%80</body>
+%72</body>
           </html>
         </notes>
       </parameter>
@@ -904,7 +834,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%81</body>
+%73</body>
           </html>
         </notes>
       </parameter>
@@ -912,7 +842,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%82</body>
+%74</body>
           </html>
         </notes>
       </parameter>
@@ -920,7 +850,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%83</body>
+%75</body>
           </html>
         </notes>
       </parameter>
@@ -928,15 +858,15 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%84</body>
+%76</body>
           </html>
         </notes>
       </parameter>
-      <parameter id="k_e1" name="k_e1" value="0" constant="true">
+      <parameter id="k_e1" name="k_e1" value="0.0005" constant="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%85</body>
+%77</body>
           </html>
         </notes>
       </parameter>
@@ -944,7 +874,7 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%86</body>
+%78</body>
           </html>
         </notes>
       </parameter>
@@ -952,50 +882,58 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%87</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="K_l0" name="K_l0" value="225.59" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%88</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="V_l0" name="V_l0" value="2.5" constant="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-%89</body>
+%79</body>
           </html>
         </notes>
       </parameter>
-      <parameter id="V_lX" name="V_lX" value="0.047" constant="true">
+      <parameter id="BW" name="BW" value="78" constant="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <body>
-%90</body>
+%80</body>
           </html>
         </notes>
       </parameter>
+      <parameter id="INS_offset" name="INS_offset" value="7" constant="true"/>
+      <parameter id="be" name="be" value="3" constant="true"/>
+      <parameter id="bradykinin" name="bradykinin" value="1" constant="true"/>
       <parameter id="pf" name="pf" value="34" constant="true"/>
-      <parameter id="U_ii" name="U_ii" value="0.8" constant="true"/>
-      <parameter id="part_m" name="part_m" value="0.5" constant="true"/>
-      <parameter id="part_f" name="part_f" value="0.5" constant="true"/>
-      <parameter id="part_l" name="part_l" value="1.7" constant="true"/>
+      <parameter id="bf_b" name="bf_b" value="3" constant="true"/>
+      <parameter id="INS_b" name="INS_b" value="0" constant="true"/>
+      <parameter id="p_bf" name="p_bf" value="1" constant="true"/>
+      <parameter id="U_ii" name="U_ii" value="0.830872" constant="true"/>
+      <parameter id="k2" name="k2" value="0.0428836" constant="true"/>
+      <parameter id="k1" name="k1" value="0.0476115" constant="true"/>
+      <parameter id="V_m" name="V_m" value="0.881325" constant="true"/>
+      <parameter id="V_mx" name="V_mx" value="0.0409359" constant="true"/>
+      <parameter id="K_m" name="K_m" value="476.351" constant="true"/>
+      <parameter id="V_l" name="V_l" value="2.00467" constant="true"/>
+      <parameter id="V_lx" name="V_lx" value="0.0439327" constant="true"/>
+      <parameter id="K_l" name="K_l" value="354.852" constant="true"/>
+      <parameter id="p1" name="p1" value="0.178911" constant="true"/>
+      <parameter id="p2" name="p2" value="4.47816" constant="true"/>
+      <parameter id="p3" name="p3" value="0.160881" constant="true"/>
+      <parameter id="p4" name="p4" value="2.62714" constant="true"/>
+      <parameter id="k_gluin" name="k_gluin" value="2.1361" constant="true"/>
+      <parameter id="k_G6P" name="k_G6P" value="11495.6" constant="true"/>
+      <parameter id="V_G6Pmax" name="V_G6Pmax" value="410.203" constant="true"/>
       <parameter id="aa" name="aa" value="0" constant="false"/>
       <parameter id="cc" name="cc" value="0" constant="false"/>
       <parameter id="EGP" name="EGP" value="0" constant="false"/>
-      <parameter id="V_mmax" name="V_mmax" value="0" constant="false"/>
-      <parameter id="V_fmax" name="V_fmax" value="0" constant="false"/>
       <parameter id="E" name="E" value="0" constant="false"/>
       <parameter id="S" name="S" value="0" constant="false"/>
       <parameter id="HE" name="HE" value="0" constant="false"/>
       <parameter id="m_3" name="m_3" value="0" constant="false"/>
       <parameter id="Ra" name="Ra" value="0" constant="false"/>
       <parameter id="k_empt" name="k_empt" value="0" constant="false"/>
+      <parameter id="bf_f" name="bf_f" value="0" constant="false"/>
+      <parameter id="bfe_f" name="bfe_f" value="0" constant="false"/>
+      <parameter id="S_po" name="S_po" value="0" constant="false"/>
+      <parameter id="V_mmax" name="V_mmax" value="0" constant="false"/>
+      <parameter id="V_lmax" name="V_lmax" value="0" constant="false"/>
+      <parameter id="INS_fe" name="INS_fe" value="0" constant="false"/>
+      <parameter id="V_in" name="V_in" value="0" constant="false"/>
+      <parameter id="V_out" name="V_out" value="0" constant="false"/>
       <parameter id="v1a" name="v1a" value="0" constant="false"/>
       <parameter id="v1basal" name="v1basal" value="0" constant="false"/>
       <parameter id="v1c" name="v1c" value="0" constant="false"/>
@@ -1030,6 +968,10 @@
       <parameter id="v9b1" name="v9b1" value="0" constant="false"/>
       <parameter id="v9f2" name="v9f2" value="0" constant="false"/>
       <parameter id="v9b2" name="v9b2" value="0" constant="false"/>
+      <parameter id="V_G6P" name="V_G6P" value="0" constant="false"/>
+      <parameter id="V_met" name="V_met" value="0" constant="false"/>
+      <parameter id="V_1" name="V_1" value="0" constant="false"/>
+      <parameter id="V_2" name="V_2" value="0" constant="false"/>
     </listOfParameters>
     <listOfRules>
       <rateRule variable="IR">
@@ -1502,6 +1444,32 @@
           </apply>
         </math>
       </rateRule>
+      <rateRule variable="Glu_in">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <apply>
+              <times/>
+              <ci> p1 </ci>
+              <apply>
+                <minus/>
+                <ci> V_in </ci>
+                <ci> V_out </ci>
+              </apply>
+            </apply>
+            <ci> V_G6P </ci>
+          </apply>
+        </math>
+      </rateRule>
+      <rateRule variable="G6P">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <ci> V_G6P </ci>
+            <ci> V_met </ci>
+          </apply>
+        </math>
+      </rateRule>
       <rateRule variable="G_p">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
@@ -1700,7 +1668,7 @@
           </apply>
         </math>
       </rateRule>
-      <rateRule variable="INS">
+      <rateRule variable="I_po">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <plus/>
@@ -1708,62 +1676,89 @@
               <times/>
               <apply>
                 <minus/>
-                <ci> p_2U </ci>
-              </apply>
-              <ci> INS </ci>
-            </apply>
-            <apply>
-              <times/>
-              <ci> p_2U </ci>
-              <apply>
-                <minus/>
-                <ci> I </ci>
-                <ci> I_b </ci>
+                <ci> gamma </ci>
               </apply>
+              <ci> I_po </ci>
             </apply>
+            <ci> S_po </ci>
           </apply>
         </math>
       </rateRule>
-      <rateRule variable="I_po">
+      <rateRule variable="Y">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
-            <plus/>
+            <times/>
             <apply>
-              <times/>
+              <minus/>
+              <ci> alpha </ci>
+            </apply>
+            <apply>
+              <minus/>
+              <ci> Y </ci>
               <apply>
-                <minus/>
-                <ci> gamma </ci>
+                <times/>
+                <ci> beta </ci>
+                <apply>
+                  <minus/>
+                  <ci> G </ci>
+                  <ci> G_b </ci>
+                </apply>
               </apply>
-              <ci> I_po </ci>
             </apply>
-            <ci> S_po </ci>
           </apply>
         </math>
       </rateRule>
-      <rateRule variable="Y">
+      <rateRule variable="INS_f">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
-            <minus/>
+            <plus/>
             <apply>
               <times/>
               <apply>
                 <minus/>
-                <ci> alpha </ci>
+                <ci> p_2U </ci>
               </apply>
-              <ci> Y </ci>
+              <ci> INS_f </ci>
             </apply>
             <apply>
               <times/>
-              <ci> beta </ci>
+              <ci> p_2U </ci>
               <apply>
                 <minus/>
-                <ci> G </ci>
-                <ci> G_b </ci>
+                <ci> I </ci>
+                <ci> I_b </ci>
               </apply>
             </apply>
           </apply>
         </math>
       </rateRule>
+      <rateRule variable="INS">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <ci> V_1 </ci>
+            <ci> V_2 </ci>
+          </apply>
+        </math>
+      </rateRule>
+      <assignmentRule variable="I">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <ci> I_p </ci>
+            <ci> V_I </ci>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="G">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <ci> G_p </ci>
+            <ci> V_G </ci>
+          </apply>
+        </math>
+      </assignmentRule>
       <assignmentRule variable="aa">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
@@ -1831,55 +1826,9 @@
           </apply>
         </math>
       </assignmentRule>
-      <assignmentRule variable="V_mmax">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <apply>
-              <minus/>
-              <cn type="integer"> 1 </cn>
-              <ci> part </ci>
-            </apply>
-            <apply>
-              <plus/>
-              <ci> V_m0 </ci>
-              <apply>
-                <times/>
-                <ci> V_mX </ci>
-                <ci> INS </ci>
-              </apply>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_fmax">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> part </ci>
-            <apply>
-              <plus/>
-              <ci> V_f0 </ci>
-              <apply>
-                <times/>
-                <ci> V_fX </ci>
-                <ci> INS </ci>
-              </apply>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
       <assignmentRule variable="E">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> k_e1 </ci>
-            <apply>
-              <minus/>
-              <ci> G_p </ci>
-              <ci> k_e2 </ci>
-            </apply>
-          </apply>
+          <cn type="integer"> 0 </cn>
         </math>
       </assignmentRule>
       <assignmentRule variable="S">
@@ -1891,24 +1840,6 @@
           </apply>
         </math>
       </assignmentRule>
-      <assignmentRule variable="I">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <ci> I_p </ci>
-            <ci> V_I </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="G">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <ci> G_p </ci>
-            <ci> V_G </ci>
-          </apply>
-        </math>
-      </assignmentRule>
       <assignmentRule variable="HE">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
@@ -2024,65 +1955,32 @@
           </apply>
         </math>
       </assignmentRule>
-      <assignmentRule variable="U_idm">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <apply>
-              <times/>
-              <ci> part_m </ci>
-              <ci> V_mmax </ci>
-              <ci> G_t </ci>
-            </apply>
-            <apply>
-              <plus/>
-              <ci> K_m0 </ci>
-              <ci> G_t </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="U_idf">
+      <assignmentRule variable="bf_f">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <times/>
-            <ci> part_f </ci>
-            <ci> nC </ci>
+            <ci> bradykinin </ci>
             <apply>
               <plus/>
-              <apply>
-                <divide/>
-                <apply>
-                  <times/>
-                  <ci> k8 </ci>
-                  <ci> GLUT4m </ci>
-                  <ci> G_t </ci>
-                </apply>
-                <ci> pf </ci>
-              </apply>
-              <apply>
-                <divide/>
-                <apply>
-                  <times/>
-                  <ci> glut1 </ci>
-                  <ci> G_t </ci>
-                </apply>
-                <ci> pf </ci>
-              </apply>
+              <ci> be </ci>
               <apply>
                 <times/>
                 <ci> kbf </ci>
                 <apply>
                   <plus/>
-                  <ci> INS </ci>
-                  <cn type="integer"> 5 </cn>
+                  <ci> INS_f </ci>
+                  <ci> INS_offset </ci>
                 </apply>
-                <cn type="e-notation"> 1 <sep/> -3 </cn>
               </apply>
             </apply>
           </apply>
         </math>
       </assignmentRule>
+      <assignmentRule variable="bfe_f">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> bf_f </ci>
+        </math>
+      </assignmentRule>
       <assignmentRule variable="S_po">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
@@ -2129,12 +2027,117 @@
           </apply>
         </math>
       </assignmentRule>
-      <assignmentRule variable="U_idl">
+      <assignmentRule variable="V_mmax">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <plus/>
+            <ci> V_m </ci>
+            <apply>
+              <times/>
+              <ci> V_mx </ci>
+              <ci> INS </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="V_lmax">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <plus/>
+            <ci> V_l </ci>
+            <apply>
+              <times/>
+              <ci> V_lx </ci>
+              <ci> INS </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="INS_fe">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <times/>
-            <ci> part_l </ci>
-            <ci> U_idm </ci>
+            <ci> nC </ci>
+            <apply>
+              <plus/>
+              <apply>
+                <divide/>
+                <apply>
+                  <times/>
+                  <ci> k8 </ci>
+                  <ci> GLUT4m </ci>
+                </apply>
+                <ci> pf </ci>
+              </apply>
+              <apply>
+                <divide/>
+                <ci> glut1 </ci>
+                <ci> pf </ci>
+              </apply>
+              <ci> bfe_f </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="V_in">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <times/>
+            <ci> p4 </ci>
+            <ci> G_t </ci>
+            <ci> INS_fe </ci>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="V_out">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <times/>
+            <ci> p3 </ci>
+            <ci> Glu_in </ci>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="U_idm">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <apply>
+              <times/>
+              <ci> V_mmax </ci>
+              <ci> G_t </ci>
+            </apply>
+            <apply>
+              <plus/>
+              <ci> K_m </ci>
+              <ci> G_t </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="U_idl">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <apply>
+              <times/>
+              <ci> V_lmax </ci>
+              <ci> G_t </ci>
+            </apply>
+            <apply>
+              <plus/>
+              <ci> K_l </ci>
+              <ci> G_t </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="U_idf">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <ci> V_in </ci>
+            <ci> V_out </ci>
           </apply>
         </math>
       </assignmentRule>
@@ -2142,8 +2145,8 @@
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <plus/>
-            <ci> U_idm </ci>
             <ci> U_idf </ci>
+            <ci> U_idm </ci>
             <ci> U_idl </ci>
           </apply>
         </math>
@@ -2156,7 +2159,7 @@
             <ci> k1a </ci>
             <apply>
               <plus/>
-              <ci> INS </ci>
+              <ci> INS_f </ci>
               <cn type="integer"> 5 </cn>
             </apply>
             <cn type="e-notation"> 1 <sep/> -3 </cn>
@@ -2536,50 +2539,66 @@
           </apply>
         </math>
       </assignmentRule>
-      <assignmentRule variable="glucoseuptake">
+      <assignmentRule variable="V_G6P">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
-            <plus/>
+            <divide/>
             <apply>
               <times/>
-              <ci> k8 </ci>
-              <ci> GLUT4m </ci>
               <apply>
                 <divide/>
                 <apply>
                   <times/>
-                  <ci> G_t </ci>
-                  <cn type="integer"> 5 </cn>
+                  <ci> V_G6Pmax </ci>
+                  <ci> Glu_in </ci>
                 </apply>
-                <cn type="integer"> 170 </cn>
-              </apply>
-            </apply>
-            <apply>
-              <times/>
-              <ci> glut1 </ci>
-              <apply>
-                <divide/>
                 <apply>
-                  <times/>
-                  <ci> G_t </ci>
-                  <cn type="integer"> 5 </cn>
+                  <plus/>
+                  <ci> k_gluin </ci>
+                  <ci> Glu_in </ci>
                 </apply>
-                <cn type="integer"> 170 </cn>
               </apply>
+              <cn type="integer"> 1 </cn>
             </apply>
             <apply>
-              <times/>
-              <ci> kbf </ci>
-              <apply>
-                <plus/>
-                <ci> INS </ci>
-                <cn type="integer"> 5 </cn>
-              </apply>
-              <cn type="e-notation"> 1 <sep/> -3 </cn>
+              <plus/>
+              <ci> k_G6P </ci>
+              <ci> G6P </ci>
+            </apply>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="V_met">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <times/>
+            <ci> p2 </ci>
+            <ci> G6P </ci>
+          </apply>
+        </math>
+      </assignmentRule>
+      <assignmentRule variable="V_1">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <times/>
+            <ci> k1 </ci>
+            <apply>
+              <minus/>
+              <ci> I </ci>
+              <ci> I_b </ci>
             </apply>
           </apply>
         </math>
       </assignmentRule>
+      <assignmentRule variable="V_2">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <times/>
+            <ci> k2 </ci>
+            <ci> INS </ci>
+          </apply>
+        </math>
+      </assignmentRule>
     </listOfRules>
   </model>
 </sbml>
diff --git a/scripts/Models/M4fat.txt b/scripts/Models/M4fat.txt
new file mode 100644
index 0000000000000000000000000000000000000000..95667ca87b4e8fc862f4823a529544c6b28b1e47
--- /dev/null
+++ b/scripts/Models/M4fat.txt
@@ -0,0 +1,213 @@
+********** MODEL NAME
+M4fat
+
+********** MODEL NOTES
+fat sub-model version of M4
+
+********** MODEL STATES
+
+d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
+d/dt(IRp)               = v1basal+v1c-v1d-v1g   %%2
+d/dt(IRins)             = v1a-v1c               %%3
+d/dt(IRip)              = v1d-v1e               %%4
+d/dt(IRi)               = v1e-v1r               %%5
+d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
+d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
+d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9 
+d/dt(X)                 = v3b-v3a               %%10
+d/dt(Xp)                = v3a-v3b               %%12
+d/dt(PKB)               = -v4a+v4b+v4h          %%13   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%14
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%15
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%16
+d/dt(mTORC1)            = v5b-v5a               %%17         
+d/dt(mTORC1a)           = v5a-v5b               %%18
+d/dt(mTORC2)            = -v5c+v5d              %%19
+d/dt(mTORC2a)           = v5c-v5d               %%20
+d/dt(AS160)             = v6b1-v6f1             %%21
+d/dt(AS160p)            = v6f1-v6b1             %%22
+d/dt(GLUT4m)            = (v7f-v7b)             %%23
+d/dt(GLUT4)             = -v7f+v7b              %%24
+d/dt(S6K)               = v9b1-v9f1             %%26
+d/dt(S6Kp)              = v9f1-v9b1             %%27
+d/dt(S6)                = v9b2-v9f2             %%28
+d/dt(S6p)               = v9f2-v9b2             %%29  
+d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
+d/dt(G6P) = V_G6P-V_met   
+ 
+IR(0)               = 100
+IRp(0)              = 0
+IRins(0)            = 0         
+IRip(0)             = 0
+IRi(0)              = 0 
+IRS1(0)             = 100
+IRS1p(0)            = 0
+IRS1p307(0)         = 0
+IRS1307(0)          = 0
+X(0)                = 100
+Xp(0)               = 0
+PKB(0)              = 100
+PKB308p(0)          = 0
+PKB473p(0)          = 0
+PKB308p473p(0)      = 0
+mTORC1(0)           = 100
+mTORC1a(0)          = 0
+mTORC2(0)           = 100
+mTORC2a(0)          = 0
+AS160(0)            = 100
+AS160p(0)           = 0
+GLUT4m(0)           = 0
+GLUT4(0)            = 100
+S6K(0)              = 100
+S6Kp(0)             = 0
+S6(0)               = 100
+S6p(0)              = 0  
+Glu_in(0) = 0
+G6P(0) = 0
+
+********** MODEL PARAMETERS
+            
+diabetes = 1  %%1                                                                
+k1a = 0.633141  %%2                                                              
+k1basal = 0.0331338  %%3                                                         
+k1c = 0.876805  %%4                                                              
+k1d = 31.012  %%5                                                                
+k1f = 1839.58  %%6                                                               
+k1g = 1944.11  %%7                                                               
+k1r = 0.547061  %%8                                                              
+k2a = 3.22728  %%9                                                               
+k2c = 5758.78  %%10                                                              
+k2basal = 0.0422768  %%11                                                        
+k2b = 3424.35  %%12                                                              
+k2d = 280.753  %%13                                                              
+k2f = 2.9131  %%14                                                               
+k2g = 0.267089  %%15                                                             
+k3a = 0.00137731  %%16                                                           
+k3b = 0.0987558  %%17                                                            
+k4a = 5790.17  %%18                                                              
+k4b = 34.7965  %%19                                                              
+k4c = 4.45581  %%20                                                              
+k4e = 42.8395  %%21                                                              
+k4f = 143.597  %%22                                                              
+k4h = 0.536145  %%23                                                             
+k5a1 = 1.8423  %%24                                                              
+k5a2 = 0.055064  %%25                                                            
+k5b = 24.826  %%26                                                               
+k5d = 1.06013  %%27                                                              
+km5 = 2.64988  %%28                                                              
+k5c = 0.0857515  %%29                                                            
+k6f1 = 2.65168  %%30                                                             
+k6f2 = 36.9348  %%31                                                             
+km6 = 30.5424  %%32                                                              
+n6 = 2.13707  %%33                                                               
+k6b = 65.1841  %%34                                                              
+k7f = 50.9829  %%35                                                              
+k7b = 2285.97  %%36                                                              
+k8 = 724.242  %%37                                                               
+glut1 = 7042.19  %%38                                                            
+k9f1 = 0.12981  %%39                                                             
+k9b1 = 0.0444092  %%40                                                           
+k9f2 = 3.3289  %%41                                                              
+k9b2 = 30.9967  %%42                                                             
+km9 = 5872.68  %%43                                                              
+n9 = 0.985466  %%44                                                              
+kbf = 0.01 %%45      1e+06                                                           
+nC = 2.1e-06  %%46  
+
+bradykinin = 1   
+INS_offset = 7 
+be = 3                                                            
+pf = 34
+
+INS1 = 1
+INS2 = 1
+INS3 = 1
+INS4 = 1
+INS5 = 1
+INS6 = 1
+INS7 = 1
+INS8 = 1
+INS9 = 1
+
+Gt1 = 1
+Gt2 = 1 
+Gt3 = 1 
+Gt4 = 1 
+Gt5 = 1 
+Gt6 = 1 
+Gt7 = 1 
+Gt8 = 1 
+Gt9 = 1
+ 
+p1 = 1 
+p2 = 2.4 
+p3 = 0.7
+p4 = 1.8 
+k_gluin = 1300
+k_G6P = 15 
+V_G6Pmax = 90 
+
+********** MODEL VARIABLES
+
+INS = interpcsIQM([0,30,60,90,120,180,240,300,360],[INS1,INS2,INS3,INS4,INS5,INS6,INS7,INS8,INS9],time)
+G_t = interpcsIQM([0,30,60,90,120,180,240,300,360],[Gt1,Gt2,Gt3,Gt4,Gt5,Gt6,Gt7,Gt8,Gt9],time)
+
+bf_f = bradykinin*(be+kbf*(INS+INS_offset))
+bfe_f = bf_f 
+
+INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
+V_in = p4*G_t*INS_fe
+V_out = p3*Glu_in
+
+U_idf = V_in-V_out
+
+
+********** MODEL REACTIONS                                                                                                                                
+
+v1a      = IR*k1a*(INS+5)*1e-3
+v1basal  = k1basal*IR
+v1c      = IRins*k1c
+v1d      = IRp*k1d
+v1e      = IRip*k1f*Xp
+v1g      = IRp*k1g
+v1r      = IRi*k1r
+v2a      = IRS1*k2a*IRip
+v2b      = IRS1p*k2b
+v2c      = IRS1p*k2c*mTORC1a*diabetes
+v2d      = IRS1p307*k2d
+v2f 	 = IRS1p307*k2f
+v2basal  = IRS1*k2basal
+v2g 	 = IRS1307*k2g
+v3a      = X*k3a*IRS1p
+v3b      = Xp*k3b
+
+v5a      = mTORC1*(k5a1*PKB308p473p+k5a2*PKB308p)
+v5b      = mTORC1a*k5b
+v5c      = mTORC2*k5c*IRip
+v5d      = k5d*mTORC2a
+v4a      = k4a*PKB*IRS1p
+v4b      = k4b*PKB308p
+v4c      = k4c*PKB308p*mTORC2a
+v4e      = k4e*PKB473p*IRS1p307
+v4f      = k4f*PKB308p473p
+v4h      = k4h*PKB473p
+
+v6f1     = AS160*(k6f1*PKB308p473p+k6f2*PKB473p^n6/(km6^n6+PKB473p^n6))
+v6b1     = AS160p*k6b
+v7f      = GLUT4*k7f*AS160p
+v7b      = GLUT4m*k7b
+
+v9f1     = S6K*k9f1*mTORC1a^n9/(km9^n9+mTORC1a^n9)
+v9b1     = S6Kp*k9b1
+v9f2     = S6*k9f2*S6Kp
+v9b2     = S6p*k9b2
+
+V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
+V_met = p2*G6P
+
+********** MODEL FUNCTIONS
+
+********** MODEL EVENTS
+
+********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/Models/M4m.txt b/scripts/Models/M4m.txt
deleted file mode 100644
index 2f7a84705df6f8811443a3cac5927e6f9ee7bbf7..0000000000000000000000000000000000000000
--- a/scripts/Models/M4m.txt
+++ /dev/null
@@ -1,127 +0,0 @@
-********** MODEL NAME
-M4m
-
-********** MODEL NOTES
-sub-model version of M4, with EGP, Ra, GLUT4m, and INS as input curves
-
-********** MODEL STATES
-          
-d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
-d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t    
-d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
-d/dt(G6P) = V_G6P-V_met   
-
-G_p(0) = 178                                                                   
-G_t(0) = 130  
-Glu_in(0) = 0
-G6P(0) = 0
-
-********** MODEL PARAMETERS
-            
-bradykinin = 1   
-kbf = 0.01
-INS_offset = 7 
-be = 3                                                            
-nC = 2.1e-06 
-k8 = 724.242                                        
-glut1 = 7042.19
-k_1 = 0.065
-k_2 = 0.079 
-pf = 34
-INS_b = 0
-p_bf = 1
-
-EGP1 = 1
-EGP2 = 1
-EGP3 = 1
-EGP4 = 1
-EGP5 = 1
-EGP6 = 1
-EGP7 = 1
-EGP8 = 1
-EGP9 = 1
-
-Ra1 = 1
-Ra2 = 1
-Ra3 = 1
-Ra4 = 1
-Ra5 = 1
-Ra6 = 1
-Ra7 = 1
-Ra8 = 1
-Ra9 = 1
-
-INS1 = 1
-INS2 = 1
-INS3 = 1
-INS4 = 1
-INS5 = 1
-INS6 = 1
-INS7 = 1
-INS8 = 1
-INS9 = 1
-
-GLm1 = 1
-GLm2 = 1 
-GLm3 = 1 
-GLm4 = 1 
-GLm5 = 1 
-GLm6 = 1 
-GLm7 = 1 
-GLm8 = 1 
-GLm9 = 1
- 
-V_m = 2.5 
-V_mx = 0.047 
-K_m = 225.59 
-part_l = 0.5
-K_l = 225.59    
-p1 = 1 
-p2 = 2.4 
-p3 = 0.7
-p4 = 1.8 
-part_f = 14  
-k_gluin = 1300
-k_G6P = 15 
-V_G6Pmax = 90 
-U_ii = 1.23
-
-********** MODEL VARIABLES
-
-INS = interpcsIQM([0,30,60,90,120,180,240,300,360],[INS1,INS2,INS3,INS4,INS5,INS6,INS7,INS8,INS9],time)
-EGP = interpcsIQM([0,30,60,90,120,180,240,300,360],[EGP1,EGP2,EGP3,EGP4,EGP5,EGP6,EGP7,EGP8,EGP9],time)
-Ra = interpcsIQM([0,30,60,90,120,180,240,300,360],[Ra1,Ra2,Ra3,Ra4,Ra5,Ra6,Ra7,Ra8,Ra9],time)
-GLUT4m = interpcsIQM([0,30,60,90,120,180,240,300,360],[GLm1,GLm2,GLm3,GLm4,GLm5,GLm6,GLm7,GLm8,GLm9],time)
-
-E = 0
-
-bf_f = bradykinin*(be+kbf*(INS+INS_offset))
-bfe_f = bf_f 
-
-part_m = 1-part_l                                     
-V_mmax = part_m*(V_m+V_mx*INS)
-V_lmax = part_l*V_mmax %(V_l+V_lx*INS) 
-
-INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
-V_in = p4*G_t*INS_fe
-V_out = p3*Glu_in
-
-U_idm = V_mmax*G_t/(K_m+G_t)
-U_idl = V_lmax*G_t/(K_l+G_t)
-U_idf = part_f*(V_in-V_out) 
-
-U_id = U_idf + U_idm + U_idl
-
-V_G = 1.88
-G = G_p/V_G
-
-********** MODEL REACTIONS                                                                                                                                
-
-V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
-V_met = p2*G6P
-
-********** MODEL FUNCTIONS
-
-********** MODEL EVENTS
-
-********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/Models/M4m.xml b/scripts/Models/M4m.xml
deleted file mode 100644
index fb22ef5c5b3c5b24245f48900b54bf4e75d330cb..0000000000000000000000000000000000000000
--- a/scripts/Models/M4m.xml
+++ /dev/null
@@ -1,438 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<sbml xmlns="http://www.sbml.org/sbml/level2" level="2" version="1">
-  <model name="M4m_SBML">
-    <notes>
-      <html xmlns="http://www.w3.org/1999/xhtml">
-        <body>
-          <p>sub-model version of M4, with EGP, Ra, GLUT4m, and INS as input curves
-</p>
-          <p>adipocyte = 1 {isCompartment:}
-</p>
-          <p>liver = 1 {isCompartment:}
-</p>
-          <p>tissue = 1 {isCompartment:}
-</p>
-          <p>muscle = 1 {isCompartment:}
-</p>
-          <p>plasma = 1 {isCompartment:}
-</p>
-          <p>fat = 1 {isCompartment:}</p>
-        </body>
-      </html>
-    </notes>
-    <listOfCompartments>
-      <compartment id="rootCompartment" name="rootCompartment" spatialDimensions="3" size="1" constant="true"/>
-    </listOfCompartments>
-    <listOfSpecies>
-      <species id="G_p" name="G_p" compartment="tissue" initialAmount="178" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="G_t" name="G_t" compartment="tissue" initialAmount="130" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="Glu_in" name="Glu_in" compartment="adipocyte" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="G6P" name="G6P" compartment="adipocyte" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="INS" name="INS" compartment="tissue" initialAmount="0" hasOnlySubstanceUnits="true" boundaryCondition="true" constant="false"/>
-      <species id="GLUT4m" name="GLUT4m" compartment="adipocyte" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_idm" name="U_idm" compartment="muscle" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_idl" name="U_idl" compartment="liver" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_idf" name="U_idf" compartment="fat" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="U_id" name="U_id" compartment="tissue" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-      <species id="G" name="G" compartment="plasma" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false"/>
-    </listOfSpecies>
-    <listOfParameters>
-      <parameter id="bradykinin" name="bradykinin" value="1" constant="true"/>
-      <parameter id="kbf" name="kbf" value="0.01" constant="true"/>
-      <parameter id="INS_offset" name="INS_offset" value="7" constant="true"/>
-      <parameter id="be" name="be" value="3" constant="true"/>
-      <parameter id="nC" name="nC" value="2.1e-006" constant="true"/>
-      <parameter id="k8" name="k8" value="724.242" constant="true"/>
-      <parameter id="glut1" name="glut1" value="7042.19" constant="true"/>
-      <parameter id="k_1" name="k_1" value="0.065" constant="true"/>
-      <parameter id="k_2" name="k_2" value="0.079" constant="true"/>
-      <parameter id="pf" name="pf" value="34" constant="true"/>
-      <parameter id="V_m" name="V_m" value="2553.38" constant="true"/>
-      <parameter id="V_mx" name="V_mx" value="2.63167" constant="true"/>
-      <parameter id="K_m" name="K_m" value="129.391" constant="true"/>
-      <parameter id="part_l" name="part_l" value="0.00249792" constant="true"/>
-      <parameter id="K_l" name="K_l" value="0.120151" constant="true"/>
-      <parameter id="p1" name="p1" value="427.534" constant="true"/>
-      <parameter id="p2" name="p2" value="423.661" constant="true"/>
-      <parameter id="p3" name="p3" value="7e-005" constant="true"/>
-      <parameter id="p4" name="p4" value="4389.24" constant="true"/>
-      <parameter id="part_f" name="part_f" value="0.932088" constant="true"/>
-      <parameter id="k_gluin" name="k_gluin" value="0.159629" constant="true"/>
-      <parameter id="k_G6P" name="k_G6P" value="0.0015" constant="true"/>
-      <parameter id="V_G6Pmax" name="V_G6Pmax" value="11.549" constant="true"/>
-      <parameter id="U_ii" name="U_ii" value="1.02548" constant="true"/>
-      <parameter id="EGP1" name="EGP1" value="1" constant="true"/>
-      <parameter id="EGP2" name="EGP2" value="1" constant="true"/>
-      <parameter id="EGP3" name="EGP3" value="1" constant="true"/>
-      <parameter id="EGP4" name="EGP4" value="1" constant="true"/>
-      <parameter id="EGP5" name="EGP5" value="1" constant="true"/>
-      <parameter id="EGP6" name="EGP6" value="1" constant="true"/>
-      <parameter id="EGP7" name="EGP7" value="1" constant="true"/>
-      <parameter id="EGP8" name="EGP8" value="1" constant="true"/>
-      <parameter id="EGP9" name="EGP9" value="1" constant="true"/>
-      <parameter id="Ra1" name="Ra1" value="1" constant="true"/>
-      <parameter id="Ra2" name="Ra2" value="1" constant="true"/>
-      <parameter id="Ra3" name="Ra3" value="1" constant="true"/>
-      <parameter id="Ra4" name="Ra4" value="1" constant="true"/>
-      <parameter id="Ra5" name="Ra5" value="1" constant="true"/>
-      <parameter id="Ra6" name="Ra6" value="1" constant="true"/>
-      <parameter id="Ra7" name="Ra7" value="1" constant="true"/>
-      <parameter id="Ra8" name="Ra8" value="1" constant="true"/>
-      <parameter id="Ra9" name="Ra9" value="1" constant="true"/>
-      <parameter id="INS1" name="INS1" value="1" constant="true"/>
-      <parameter id="INS2" name="INS2" value="1" constant="true"/>
-      <parameter id="INS3" name="INS3" value="1" constant="true"/>
-      <parameter id="INS4" name="INS4" value="1" constant="true"/>
-      <parameter id="INS5" name="INS5" value="1" constant="true"/>
-      <parameter id="INS6" name="INS6" value="1" constant="true"/>
-      <parameter id="INS7" name="INS7" value="1" constant="true"/>
-      <parameter id="INS8" name="INS8" value="1" constant="true"/>
-      <parameter id="INS9" name="INS9" value="1" constant="true"/>
-      <parameter id="GLm1" name="GLm1" value="1" constant="true"/>
-      <parameter id="GLm2" name="GLm2" value="1" constant="true"/>
-      <parameter id="GLm3" name="GLm3" value="1" constant="true"/>
-      <parameter id="GLm4" name="GLm4" value="1" constant="true"/>
-      <parameter id="GLm5" name="GLm5" value="1" constant="true"/>
-      <parameter id="GLm6" name="GLm6" value="1" constant="true"/>
-      <parameter id="GLm7" name="GLm7" value="1" constant="true"/>
-      <parameter id="GLm8" name="GLm8" value="1" constant="true"/>
-      <parameter id="GLm9" name="GLm9" value="1" constant="true"/>
-      <parameter id="EGP" name="EGP" value="0" constant="false"/>
-      <parameter id="Ra" name="Ra" value="0" constant="false"/>
-      <parameter id="E" name="E" value="0" constant="false"/>
-      <parameter id="bf_f" name="bf_f" value="0" constant="false"/>
-      <parameter id="bfe_f" name="bfe_f" value="0" constant="false"/>
-      <parameter id="part_m" name="part_m" value="0" constant="false"/>
-      <parameter id="V_mmax" name="V_mmax" value="0" constant="false"/>
-      <parameter id="V_lmax" name="V_lmax" value="0" constant="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-(V_l+V_lx*INS)</body>
-          </html>
-        </notes>
-      </parameter>
-      <parameter id="INS_fe" name="INS_fe" value="0" constant="false"/>
-      <parameter id="V_in" name="V_in" value="0" constant="false"/>
-      <parameter id="V_out" name="V_out" value="0" constant="false"/>
-      <parameter id="V_G" name="V_G" value="0" constant="false"/>
-      <parameter id="V_G6P" name="V_G6P" value="0" constant="false"/>
-      <parameter id="V_met" name="V_met" value="0" constant="false"/>
-    </listOfParameters>
-    <listOfRules>
-      <rateRule variable="G_p">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <plus/>
-            <apply>
-              <minus/>
-              <apply>
-                <minus/>
-                <apply>
-                  <minus/>
-                  <apply>
-                    <plus/>
-                    <ci> EGP </ci>
-                    <ci> Ra </ci>
-                  </apply>
-                  <ci> E </ci>
-                </apply>
-                <ci> U_ii </ci>
-              </apply>
-              <apply>
-                <times/>
-                <ci> k_1 </ci>
-                <ci> G_p </ci>
-              </apply>
-            </apply>
-            <apply>
-              <times/>
-              <ci> k_2 </ci>
-              <ci> G_t </ci>
-            </apply>
-          </apply>
-        </math>
-      </rateRule>
-      <rateRule variable="G_t">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <minus/>
-            <apply>
-              <plus/>
-              <apply>
-                <minus/>
-                <ci> U_id </ci>
-              </apply>
-              <apply>
-                <times/>
-                <ci> k_1 </ci>
-                <ci> G_p </ci>
-              </apply>
-            </apply>
-            <apply>
-              <times/>
-              <ci> k_2 </ci>
-              <ci> G_t </ci>
-            </apply>
-          </apply>
-        </math>
-      </rateRule>
-      <rateRule variable="Glu_in">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <minus/>
-            <apply>
-              <times/>
-              <ci> p1 </ci>
-              <apply>
-                <minus/>
-                <ci> V_in </ci>
-                <ci> V_out </ci>
-              </apply>
-            </apply>
-            <ci> V_G6P </ci>
-          </apply>
-        </math>
-      </rateRule>
-      <rateRule variable="G6P">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <minus/>
-            <ci> V_G6P </ci>
-            <ci> V_met </ci>
-          </apply>
-        </math>
-      </rateRule>
-      <assignmentRule variable="INS"/>
-      <assignmentRule variable="EGP"/>
-      <assignmentRule variable="Ra"/>
-      <assignmentRule variable="GLUT4m"/>
-      <assignmentRule variable="E">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <cn type="integer"> 0 </cn>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="bf_f">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> bradykinin </ci>
-            <apply>
-              <plus/>
-              <ci> be </ci>
-              <apply>
-                <times/>
-                <ci> kbf </ci>
-                <apply>
-                  <plus/>
-                  <ci> INS </ci>
-                  <ci> INS_offset </ci>
-                </apply>
-              </apply>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="bfe_f">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <ci> bf_f </ci>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="part_m">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <minus/>
-            <cn type="integer"> 1 </cn>
-            <ci> part_l </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_mmax">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> part_m </ci>
-            <apply>
-              <plus/>
-              <ci> V_m </ci>
-              <apply>
-                <times/>
-                <ci> V_mx </ci>
-                <ci> INS </ci>
-              </apply>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_lmax">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <body>
-(V_l+V_lx*INS)</body>
-          </html>
-        </notes>
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> part_l </ci>
-            <ci> V_mmax </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="INS_fe">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> nC </ci>
-            <apply>
-              <plus/>
-              <apply>
-                <divide/>
-                <apply>
-                  <times/>
-                  <ci> k8 </ci>
-                  <ci> GLUT4m </ci>
-                </apply>
-                <ci> pf </ci>
-              </apply>
-              <apply>
-                <divide/>
-                <ci> glut1 </ci>
-                <ci> pf </ci>
-              </apply>
-              <ci> bfe_f </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_in">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> p4 </ci>
-            <ci> G_t </ci>
-            <ci> INS_fe </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_out">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> p3 </ci>
-            <ci> Glu_in </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="U_idm">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <apply>
-              <times/>
-              <ci> V_mmax </ci>
-              <ci> G_t </ci>
-            </apply>
-            <apply>
-              <plus/>
-              <ci> K_m </ci>
-              <ci> G_t </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="U_idl">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <apply>
-              <times/>
-              <ci> V_lmax </ci>
-              <ci> G_t </ci>
-            </apply>
-            <apply>
-              <plus/>
-              <ci> K_l </ci>
-              <ci> G_t </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="U_idf">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> part_f </ci>
-            <apply>
-              <minus/>
-              <ci> V_in </ci>
-              <ci> V_out </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="U_id">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <plus/>
-            <ci> U_idf </ci>
-            <ci> U_idm </ci>
-            <ci> U_idl </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_G">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <cn> 1.88 </cn>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="G">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <ci> G_p </ci>
-            <ci> V_G </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_G6P">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <divide/>
-            <apply>
-              <times/>
-              <apply>
-                <divide/>
-                <apply>
-                  <times/>
-                  <ci> V_G6Pmax </ci>
-                  <ci> Glu_in </ci>
-                </apply>
-                <apply>
-                  <plus/>
-                  <ci> k_gluin </ci>
-                  <ci> Glu_in </ci>
-                </apply>
-              </apply>
-              <cn type="integer"> 1 </cn>
-            </apply>
-            <apply>
-              <plus/>
-              <ci> k_G6P </ci>
-              <ci> G6P </ci>
-            </apply>
-          </apply>
-        </math>
-      </assignmentRule>
-      <assignmentRule variable="V_met">
-        <math xmlns="http://www.w3.org/1998/Math/MathML">
-          <apply>
-            <times/>
-            <ci> p2 </ci>
-            <ci> G6P </ci>
-          </apply>
-        </math>
-      </assignmentRule>
-    </listOfRules>
-  </model>
-</sbml>
diff --git a/scripts/Models/M4pred.txt b/scripts/Models/M4pred.txt
new file mode 100644
index 0000000000000000000000000000000000000000..42dfe787d4ae602f30afb7e7776e5d3d4ac7e8d0
--- /dev/null
+++ b/scripts/Models/M4pred.txt
@@ -0,0 +1,288 @@
+********** MODEL NAME
+M4pred
+********** MODEL NOTES
+M4 with initial values of the validation data
+
+********** MODEL STATES
+d/dt(IR)                = -v1a-v1basal+v1r+v1g  %%1 
+d/dt(IRp)               = v1basal+v1c-v1d-v1g   %%2
+d/dt(IRins)             = v1a-v1c               %%3
+d/dt(IRip)              = v1d-v1e               %%4
+d/dt(IRi)               = v1e-v1r               %%5
+d/dt(IRS1)              = v2b+v2g-v2a-v2basal   %%6
+d/dt(IRS1p)             = v2a+v2d-v2b-v2c       %%7
+d/dt(IRS1p307)          = v2c-v2d-v2f           %%8
+d/dt(IRS1307)           = v2basal+v2f-v2g       %%9
+d/dt(X)                 = v3b-v3a               %%10
+d/dt(Xp)                = v3a-v3b               %%11
+d/dt(PKB)               = -v4a+v4b+v4h          %%12   
+d/dt(PKB308p)           = v4a-v4b-v4c           %%13
+d/dt(PKB473p)           = -v4e+v4f-v4h          %%14
+d/dt(PKB308p473p)       = v4c+v4e-v4f           %%15
+d/dt(mTORC1)            = v5b-v5a               %%16         
+d/dt(mTORC1a)           = v5a-v5b               %%17
+d/dt(mTORC2)            = -v5c+v5d              %%18
+d/dt(mTORC2a)           = v5c-v5d               %%19
+d/dt(AS160)             = v6b1-v6f1             %%20
+d/dt(AS160p)            = v6f1-v6b1             %%21
+d/dt(GLUT4m)            = v7f-v7b               %%22
+d/dt(GLUT4)             = -v7f+v7b              %%23
+d/dt(S6K)               = v9b1-v9f1             %%24
+d/dt(S6Kp)              = v9f1-v9b1             %%25
+d/dt(S6)                = v9b2-v9f2             %%26
+d/dt(S6p)               = v9f2-v9b2             %%27             
+
+d/dt(Glu_in) = p1*(V_in-V_out)-V_G6P 
+d/dt(G6P) = V_G6P-V_met           
+                                                      
+d/dt(G_p) = EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t
+d/dt(G_t) = (-U_id)+k_1*G_p-k_2*G_t                                      
+d/dt(I_l) = (-m_1*I_l)-m_3*I_l+m_2*I_p+S  
+d/dt(I_p) = (-m_2*I_p)-m_4*I_p+m_1*I_l                                            
+d/dt(Q_sto1) = -k_gri*Q_sto1                                                     
+d/dt(Q_sto2) = (-k_empt*Q_sto2)+k_gri*Q_sto1                                     
+d/dt(Q_gut) = (-k_abs*Q_gut)+k_empt*Q_sto2                                       
+d/dt(I_1) = -k_i*(I_1-I)                                                         
+d/dt(I_d) = -k_i*(I_d-I_1)                                                           
+d/dt(I_po) = (-gamma*I_po)+S_po                                                  
+d/dt(Y) = -alpha*(Y-beta*(G-G_b))                                                        
+d/dt(INS_f) = (-p_2U*INS_f)+p_2U*(I-I_b)  
+d/dt(INS) = V_1-V_2    
+
+IR(0)               = 100
+IRp(0)              = 0
+IRins(0)            = 0         
+IRip(0)             = 0
+IRi(0)              = 0 
+IRS1(0)             = 100
+IRS1p(0)            = 0
+IRS1p307(0)         = 0
+IRS1307(0)          = 0
+X(0)                = 100
+Xp(0)               = 0
+PKB(0)              = 100
+PKB308p(0)          = 0
+PKB473p(0)          = 0
+PKB308p473p(0)      = 0
+mTORC1(0)           = 100
+mTORC1a(0)          = 0
+mTORC2(0)           = 100
+mTORC2a(0)          = 0
+AS160(0)            = 100
+AS160p(0)           = 0
+GLUT4m(0)           = 0
+GLUT4(0)            = 100
+S6K(0)              = 100
+S6Kp(0)             = 0
+S6(0)               = 100
+S6p(0)              = 0    
+
+Glu_in(0) = 0
+G6P(0) = 0
+                           
+G_p(0) = 167.4189                                                                   
+G_t(0) = 130                                                                     
+I_l(0) = 4.5                                                                 
+I_p(0) = 1.5873                                                                 
+Q_sto1(0) = 78000                                                                
+Q_sto2(0) = 0                                                                    
+Q_gut(0) = 0                                                                     
+I_1(0) = 25                                                                      
+I_d(0) = 25                                                                      
+I_po(0) = 3.6                                                  
+Y(0) = 0
+INS_f(0) = 0
+INS(0) = 0
+
+********** MODEL PARAMETERS 
+
+diabetes = 1  %%1                                                                
+k1a = 0.633141  %%2                                                              
+k1basal = 0.0331338  %%3                                                         
+k1c = 0.876805  %%4                                                              
+k1d = 31.012  %%5                                                                
+k1f = 1839.58  %%6                                                               
+k1g = 1944.11  %%7                                                               
+k1r = 0.547061  %%8                                                              
+k2a = 3.22728  %%9                                                               
+k2c = 5758.78  %%10                                                              
+k2basal = 0.0422768  %%11                                                        
+k2b = 3424.35  %%12                                                              
+k2d = 280.753  %%13                                                              
+k2f = 2.9131  %%14                                                               
+k2g = 0.267089  %%15                                                             
+k3a = 0.00137731  %%16                                                           
+k3b = 0.0987558  %%17                                                            
+k4a = 5790.17  %%18                                                              
+k4b = 34.7965  %%19                                                              
+k4c = 4.45581  %%20                                                              
+k4e = 42.8395  %%21                                                              
+k4f = 143.597  %%22                                                              
+k4h = 0.536145  %%23                                                             
+k5a1 = 1.8423  %%24                                                              
+k5a2 = 0.055064  %%25                                                            
+k5b = 24.826  %%26                                                               
+k5d = 1.06013  %%27                                                              
+km5 = 2.64988  %%28                                                              
+k5c = 0.0857515  %%29                                                            
+k6f1 = 2.65168  %%30                                                             
+k6f2 = 36.9348  %%31                                                             
+km6 = 30.5424  %%32                                                              
+n6 = 2.13707  %%33                                                               
+k6b = 65.1841  %%34                                                              
+k7f = 50.9829  %%35                                                              
+k7b = 2285.97  %%36                                                              
+k8 = 724.242  %%37                                                               
+glut1 = 7042.19  %%38                                                            
+k9f1 = 0.12981  %%39                                                             
+k9b1 = 0.0444092  %%40                                                           
+k9f2 = 3.3289  %%41                                                              
+k9b2 = 30.9967  %%42                                                             
+km9 = 5872.68  %%43                                                              
+n9 = 0.985466  %%44                                                              
+kbf = 0.01 %%45      1e+06                                                           
+nC = 2.1e-06  %%46  
+
+V_G = 1.88  %%47                                                                 
+k_1 = 0.065  %%48                                                                
+k_2 = 0.079  %%49                                                                
+G_b = 95  %%50                                                                   
+V_I = 0.05  %%51                                                                 
+m_1 = 0.19  %%52                                                                 
+m_2 = 0.484    %% 53                                                                  
+m_4 = 0.194    %% 54                                                                  
+m_5 = 0.0304  %%55                                                               
+m_6 = 0.6471   %% 56                                                                  
+HE_b = 0.6  %%57                                                                 
+I_b = 25  %%58                                                                   
+S_b = 1.8  %%59                                                                  
+k_max = 0.0558  %%60                                                             
+k_min = 0.008  %%61                                                              
+k_abs = 0.057  %%62                                                              
+k_gri = 0.0558  %%63                                                             
+f = 0.9  %%64                                                                    
+b = 0.82  %%65                                                                   
+d = 0.01  %%66                                                                   
+k_p1 = 2.7  %%68                                                                 
+k_p2 = 0.0021  %%69                                                              
+k_p3 = 0.009  %%70                                                               
+k_p4 = 0.0618  %%71                                                              
+k_i = 0.0079  %%72                                                               
+p_2U = 0.0331  %%80                                                              
+part = 0.2  %%81                                                                 
+K = 2.3  %%82                                                                    
+alpha = 0.05  %%83                                                               
+beta = 0.11  %%84                                                                
+gamma = 0.5  %%85                                                                
+k_e1 = 0.0005  %%86                                                              
+k_e2 = 339  %%87
+D = 75000  %%88                                                                   
+BW = 75    %%89
+
+INS_offset = 7 
+be = 3          
+bradykinin = 1  
+pf = 34
+bf_b = 3 
+INS_b = 0                                                               
+p_bf = 1                                                     
+                                            
+U_ii = 1.47  
+k2 = 0.0331 
+k1 = 0.0331
+V_m = 2.5 
+V_mx = 0.047 
+K_m = 225.59 
+V_l = 2.5 
+V_lx = 0.047 
+K_l = 225.59
+p1 = 1 
+p2 = 2.4 
+p3 = 0.7
+p4 = 1.8  
+k_gluin = 1300
+k_G6P = 15 
+V_G6Pmax = 90  
+
+********** MODEL VARIABLES
+
+I = I_p/V_I                                                                      
+G = G_p/V_G  
+aa = 5/2/(1-b)/D                                                                 
+cc = 5/2/d/D                                                                     
+EGP = k_p1-k_p2*G_p-k_p3*I_d-k_p4*I_po                                                      
+E = 0                                                                            
+S = gamma*I_po                                                                     
+HE = (-m_5*S)+m_6                                                                
+m_3 = HE*m_1/(1-HE)                                                              
+Q_sto = Q_sto1+Q_sto2                                                            
+Ra = f*k_abs*Q_gut/BW                                                           
+k_empt = k_min+(k_max-k_min)/2*(tanh(aa*(Q_sto-b*D))-tanh(cc*(Q_sto-d*D))+2)   
+bf_f = bradykinin*(be+kbf*(INS_f+INS_offset))
+bfe_f = bf_f %(bf_f-bf_b)*(INS_f-INS_b)*p_bf
+S_po = Y+K*(EGP+Ra-E-U_ii-k_1*G_p+k_2*G_t)/V_G+S_b
+                                  
+V_mmax = V_m+V_mx*INS
+V_lmax = V_l+V_lx*INS
+
+INS_fe = nC*(k8*GLUT4m/pf + glut1/pf + bfe_f)
+V_in = p4*G_t*INS_fe
+V_out = p3*Glu_in
+
+U_idm = V_mmax*G_t/(K_m+G_t)
+U_idl = V_lmax*G_t/(K_l+G_t)
+U_idf = V_in-V_out 
+
+U_id = U_idf + U_idm + U_idl
+
+********** MODEL REACTIONS                                                                                                                                
+
+v1a      = IR*k1a*(INS_f+5)*1e-3
+v1basal  = k1basal*IR
+v1c      = IRins*k1c
+v1d      = IRp*k1d
+v1e      = IRip*k1f*Xp
+v1g      = IRp*k1g
+v1r      = IRi*k1r
+v2a      = IRS1*k2a*IRip
+v2b      = IRS1p*k2b
+v2c      = IRS1p*k2c*mTORC1a*diabetes
+v2d      = IRS1p307*k2d
+v2f 	 = IRS1p307*k2f
+v2basal  = IRS1*k2basal
+v2g 	 = IRS1307*k2g
+v3a      = X*k3a*IRS1p
+v3b      = Xp*k3b
+
+v5a      = mTORC1*(k5a1*PKB308p473p+k5a2*PKB308p)
+v5b      = mTORC1a*k5b
+v5c      = mTORC2*k5c*IRip
+v5d      = k5d*mTORC2a
+v4a      = k4a*PKB*IRS1p
+v4b      = k4b*PKB308p
+v4c      = k4c*PKB308p*mTORC2a
+v4e      = k4e*PKB473p*IRS1p307
+v4f      = k4f*PKB308p473p
+v4h      = k4h*PKB473p
+
+v6f1     = AS160*(k6f1*PKB308p473p+k6f2*PKB473p^n6/(km6^n6+PKB473p^n6))
+v6b1     = AS160p*k6b
+v7f      = GLUT4*k7f*AS160p
+v7b      = GLUT4m*k7b
+
+v9f1     = S6K*k9f1*mTORC1a^n9/(km9^n9+mTORC1a^n9)
+v9b1     = S6Kp*k9b1
+v9f2     = S6*k9f2*S6Kp
+v9b2     = S6p*k9b2
+
+V_G6P = V_G6Pmax*Glu_in/(k_gluin+Glu_in)*1/(k_G6P+G6P)
+V_met = p2*G6P
+
+V_1 = k1*(I-I_b)
+V_2 = k2*INS 
+
+********** MODEL FUNCTIONS
+
+********** MODEL EVENTS
+
+********** MODEL MATLAB FUNCTIONS
\ No newline at end of file
diff --git a/scripts/costFun.m b/scripts/costFunD.m
similarity index 52%
rename from scripts/costFun.m
rename to scripts/costFunD.m
index 0cf965e0f698774db94c8b5226496b20e623e84f..9d3b5d5be2c2a35c557df565c16c33bec6a6a0af 100644
--- a/scripts/costFun.m
+++ b/scripts/costFunD.m
@@ -1,5 +1,18 @@
 
-function cost = costFun(params,DATA,model,fid,limit,paramsFixed)
+function cost = costFunD(DATA,model,params,fid,limit,paramsFixed)
+
+global SEM_EGP
+global EGP_dalla
+global SEM_S
+global S_dalla
+global SEM_I
+global pi_dalla
+global SEM_G
+global pg_dalla
+global SEM_Ra
+global Ra_dalla
+global SEM_U
+global GU_dalla
 
 try
 [sim_muscle, sim_adipose, sim_liver,sim_ii]=simulateModel(DATA.time,model,exp(params),paramsFixed);
@@ -13,11 +26,11 @@ if trapz(DATA.time,sim_liver)<0.38*(trapz(DATA.time,sim_adipose)+trapz(DATA.time
 cost=1000+abs(1/(trapz(DATA.time,sim_liver)+1));
 end
 
-if trapz(DATA.time,sim_liver)>0.55*(trapz(DATA.time,sim_adipose)+trapz(DATA.time,sim_muscle)+trapz(DATA.time,sim_ii)+trapz(DATA.time,sim_liver))
+if trapz(DATA.time,sim_liver)>0.49*(trapz(DATA.time,sim_adipose)+trapz(DATA.time,sim_muscle)+trapz(DATA.time,sim_ii)+trapz(DATA.time,sim_liver))
 cost=1000+abs(0.1*(trapz(DATA.time,sim_liver)));
 end
 
-if trapz(DATA.time,sim_ii)<0.16*(trapz(DATA.time,sim_adipose)+trapz(DATA.time,sim_muscle)+trapz(DATA.time,sim_ii)+trapz(DATA.time,sim_liver))
+if trapz(DATA.time,sim_ii)<0.18*(trapz(DATA.time,sim_adipose)+trapz(DATA.time,sim_muscle)+trapz(DATA.time,sim_ii)+trapz(DATA.time,sim_liver))
 cost=1000+abs(1/(trapz(DATA.time,sim_ii)+1));
 end
 
@@ -25,6 +38,22 @@ if trapz(DATA.time,sim_ii)>0.30*(trapz(DATA.time,sim_adipose)+trapz(DATA.time,si
 cost=1000+abs(0.1*(trapz(DATA.time,sim_ii)));
 end
 
+%%% dallaman
+timeDalla = pg_dalla(:,1)';
+[Gp,Ip,EGP,Ra,U_id,S] = simulateModelDalla(timeDalla,model,exp(params),paramsFixed);
+
+costGp = nansum((((pg_dalla(:,2)'-Gp).^2)./(SEM_G'.^2))) ;
+costIp = nansum((((pi_dalla(:,2)'-Ip).^2)./(SEM_I'.^2))) ;
+costEGP = nansum((((EGP_dalla(:,2)'-EGP).^2)./(SEM_EGP'.^2))) ;
+costRa = nansum((((Ra_dalla(:,2)'-Ra).^2)./(SEM_Ra'.^2))) ;
+U = U_id + params(1).*ones(size(Gp));
+costU = nansum((((GU_dalla(:,2)'-U).^2)./(SEM_U'.^2))) ;
+costS = nansum((((S_dalla(:,2)'-S).^2)./(SEM_S'.^2))) ;
+
+costDalla = costGp+costIp+costEGP+costRa+costU+costS;
+
+cost = cost+costDalla;
+
 if cost<limit 
  fprintf(fid, [sprintf('%.52f, ',exp(params)) sprintf('%.52f\n',cost)]); 
 end
diff --git a/scripts/costFunDVal.m b/scripts/costFunDVal.m
new file mode 100644
index 0000000000000000000000000000000000000000..9b0f3b5f3d5058e2b087a5a9f4e7bc3c97fc200b
--- /dev/null
+++ b/scripts/costFunDVal.m
@@ -0,0 +1,11 @@
+function cost = costFunDVal(model,params,paramsFixed,OGTT)
+
+
+[Gp,Ip,~,~,~,~] = simulateModelDalla(OGTT.time,model,params,paramsFixed);
+
+cost1 = nansum((((OGTT.Glucose(2:end)-Gp(2:end-1)').^2)./(OGTT.SEMg(2:end).^2)));
+cost2 = nansum((((OGTT.Insulin(2:end)-Ip(2:end)').^2)./(OGTT.SEMi(2:end).^2)));
+
+cost = cost1 + cost2;
+
+end
\ No newline at end of file
diff --git a/scripts/costFunDalla.m b/scripts/costFunDalla.m
new file mode 100644
index 0000000000000000000000000000000000000000..901691c100387b1a8a2a35b179f2a2f20f4df5be
--- /dev/null
+++ b/scripts/costFunDalla.m
@@ -0,0 +1,37 @@
+
+function cost = costFunDalla(model,params,paramsFixed)
+
+global SEM_EGP
+global EGP_dalla
+global SEM_S
+global S_dalla
+global SEM_I
+global pi_dalla
+global SEM_G
+global pg_dalla
+global SEM_Ra
+global Ra_dalla
+global SEM_U
+global GU_dalla
+
+try
+
+timeDalla = pg_dalla(:,1)';
+[Gp,Ip,EGP,Ra,U_id,S] = simulateModelDalla(timeDalla,model,params,paramsFixed);
+
+costGp = nansum((((pg_dalla(:,2)'-Gp).^2)./(SEM_G'.^2))) ;
+costIp = nansum((((pi_dalla(:,2)'-Ip).^2)./(SEM_I'.^2))) ;
+costEGP = nansum((((EGP_dalla(:,2)'-EGP).^2)./(SEM_EGP'.^2))) ;
+costRa = nansum((((Ra_dalla(:,2)'-Ra).^2)./(SEM_Ra'.^2))) ;
+U = U_id + params(1).*ones(size(Gp));
+costU = nansum((((GU_dalla(:,2)'-U).^2)./(SEM_U'.^2))) ;
+costS = nansum((((S_dalla(:,2)'-S).^2)./(SEM_S'.^2))) ;
+
+cost = costGp+costIp+costEGP+costRa+costU+costS;
+
+catch err
+disp(getReport(err))
+    cost=1e99;
+end
+
+end
diff --git a/scripts/costFunFat.m b/scripts/costFunFat.m
new file mode 100644
index 0000000000000000000000000000000000000000..1e70cc7c3e8142287d76a6a2fcb10d849a5bd9bc
--- /dev/null
+++ b/scripts/costFunFat.m
@@ -0,0 +1,18 @@
+
+function cost = costFunFat(DATA,model,params,fid,limit,paramsFixed)
+
+try
+[~,sim_adipose,~,~]=simulateModel(DATA.time,model,exp(params),paramsFixed);
+
+cost = nansum((((DATA.Adipose'-sim_adipose).^2)./(DATA.SEM2'.^2)));
+
+if cost<limit 
+ fprintf(fid, [sprintf('%.52f, ',exp(params)) sprintf('%.52f\n',cost)]); 
+end
+
+catch err
+disp(getReport(err))
+    cost=1e99;
+end
+
+end
diff --git a/scripts/costFunVal.m b/scripts/costFunVal.m
new file mode 100644
index 0000000000000000000000000000000000000000..c994553b97117d123fb2a3c06373c8283b37f640
--- /dev/null
+++ b/scripts/costFunVal.m
@@ -0,0 +1,11 @@
+
+function cost = costFunVal(DATA,model,params,paramsFixed)
+
+[sim_muscle, sim_adipose,~,~]=simulateModel(DATA.time,model,params,paramsFixed);
+ 
+cost1 = nansum((((DATA.Muscle-sim_muscle').^2)./(DATA.SEM1.^2)));
+cost2 = nansum((((DATA.Adipose-sim_adipose').^2)./(DATA.SEM2.^2)));
+
+cost = cost1 + cost2;
+
+end
diff --git a/scripts/data/AdamsHigh.mat b/scripts/data/AdamsHigh.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2d1476350f3f6b0d0b0d88c5ba33b0e873032cd6
Binary files /dev/null and b/scripts/data/AdamsHigh.mat differ
diff --git a/scripts/data/Coppack_fusk.mat b/scripts/data/Coppack_fusk.mat
new file mode 100644
index 0000000000000000000000000000000000000000..de79c4f4b4e6f10ac5f29e3351bc563162de4996
Binary files /dev/null and b/scripts/data/Coppack_fusk.mat differ
diff --git a/scripts/data/Dimitriadis.mat b/scripts/data/Dimitriadis.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b4dff515d9c73af4725a132047da1bd12b5cc139
Binary files /dev/null and b/scripts/data/Dimitriadis.mat differ
diff --git a/scripts/data/EGP.mat b/scripts/data/EGP.mat
index bfc43bd234fb622d2e6907a832a52e94133679e3..e940907bb21f93091c088bbc4a947a6d97274888 100644
Binary files a/scripts/data/EGP.mat and b/scripts/data/EGP.mat differ
diff --git a/scripts/data/GLUT4m.mat b/scripts/data/GLUT4m.mat
index 2de30ec9cfeda248cbec29ff69e70a8db75ff3dd..24f6f6d17758c1a7468db169facd59c5ef3d9140 100644
Binary files a/scripts/data/GLUT4m.mat and b/scripts/data/GLUT4m.mat differ
diff --git a/scripts/data/G_p.mat b/scripts/data/G_p.mat
new file mode 100644
index 0000000000000000000000000000000000000000..fd64117c83f040d20ff356d14ec7e4ef0b76a1ba
Binary files /dev/null and b/scripts/data/G_p.mat differ
diff --git a/scripts/data/G_t.mat b/scripts/data/G_t.mat
index 47f6a16867aa93f153f37f4f215814b526dbe617..8d819ced934fe38f541029cf3fbe0dae14843b94 100644
Binary files a/scripts/data/G_t.mat and b/scripts/data/G_t.mat differ
diff --git a/scripts/data/I.mat b/scripts/data/I.mat
index 7af915968668c644beb24ad61df8b75a3ca744da..e5aa193b024207a6079b251521449c5191f69ea7 100644
Binary files a/scripts/data/I.mat and b/scripts/data/I.mat differ
diff --git a/scripts/data/INS.mat b/scripts/data/INS.mat
index ed9e22e5e8ca8dd698092012c71a9f3f66154de2..0c1bca886bd57e3b70bf7d512df70cf3dc5d1c46 100644
Binary files a/scripts/data/INS.mat and b/scripts/data/INS.mat differ
diff --git a/scripts/data/Ra.mat b/scripts/data/Ra.mat
index 2d7cd7dfaf054d20a0a123ddbce55e52d62327e4..77ec3009c2df4a1a9d81be28227a707423c6a37b 100644
Binary files a/scripts/data/Ra.mat and b/scripts/data/Ra.mat differ
diff --git a/scripts/data/SEM_EGP.mat b/scripts/data/SEM_EGP.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3a1bd446f100a847a24b6487c0253009c3f7c849
Binary files /dev/null and b/scripts/data/SEM_EGP.mat differ
diff --git a/scripts/data/SEM_G.mat b/scripts/data/SEM_G.mat
new file mode 100644
index 0000000000000000000000000000000000000000..23dbe5a3a6caac6cc4b204da5af9c4369c5498a6
Binary files /dev/null and b/scripts/data/SEM_G.mat differ
diff --git a/scripts/data/SEM_I.mat b/scripts/data/SEM_I.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5557a6eb6caa8253250be106ff8e573eb302b8b8
Binary files /dev/null and b/scripts/data/SEM_I.mat differ
diff --git a/scripts/data/SEM_Ra.mat b/scripts/data/SEM_Ra.mat
new file mode 100644
index 0000000000000000000000000000000000000000..fe8c255694391e9d233106c921e4496e099fa179
Binary files /dev/null and b/scripts/data/SEM_Ra.mat differ
diff --git a/scripts/data/SEM_S.mat b/scripts/data/SEM_S.mat
new file mode 100644
index 0000000000000000000000000000000000000000..42656cd82099d3ae451115b945075cd3ff48464b
Binary files /dev/null and b/scripts/data/SEM_S.mat differ
diff --git a/scripts/data/SEM_U.mat b/scripts/data/SEM_U.mat
new file mode 100644
index 0000000000000000000000000000000000000000..86e1b7b0aae58d4cac1eec0fdd29e315c6e4a033
Binary files /dev/null and b/scripts/data/SEM_U.mat differ
diff --git a/scripts/data/U_idf.mat b/scripts/data/U_idf.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8fb6a8d0b056d0a3ffeb1026a9780a703f458d9f
Binary files /dev/null and b/scripts/data/U_idf.mat differ
diff --git a/scripts/data/pg_dalla.mat b/scripts/data/pg_dalla.mat
index 19c556cc14113d1e47b3b951a0320c33b58327b5..d58fe751c2e60c77fc7ab049b8f4b4191ad3241e 100644
Binary files a/scripts/data/pg_dalla.mat and b/scripts/data/pg_dalla.mat differ
diff --git a/scripts/data/pg_high.mat b/scripts/data/pg_high.mat
index 9a80bfd27085b38c22a142a6e8af99a585dcd966..640b1ba7ea7c1692a22752da339b255be1dcd7b2 100644
Binary files a/scripts/data/pg_high.mat and b/scripts/data/pg_high.mat differ
diff --git a/scripts/data/pg_low.mat b/scripts/data/pg_low.mat
index 7a34962a0ee38b0a4083ebb0f635e95defd6d546..bdbeee454b633236f51c731b60d95de6d12723ed 100644
Binary files a/scripts/data/pg_low.mat and b/scripts/data/pg_low.mat differ
diff --git a/scripts/getMaxMin.m b/scripts/getMaxMin.m
index 87e19be947d986d4b8caaab5bedc080846266c97..c8b0ec6c0e518beee83ea8b2b92c499b628db40c 100644
--- a/scripts/getMaxMin.m
+++ b/scripts/getMaxMin.m
@@ -1,4 +1,4 @@
-function [boundries,optParams]=getMaxMin(simTime, model, params,paramsFixed)
+function [boundries]=getMaxMin(simTime, model, params,paramsFixed,DATA,limit)
 simVal=inf(2,length(simTime));
 
 minVal=simVal;
@@ -17,13 +17,11 @@ optParams = nan(size(params(1,:)));
                 minVal=simVal;
                 maxVal=minVal;
             end
-            if max(abs(simVal(:,:)))<inf
+            if max(abs(simVal(:,:)))<inf  
+                if costFunVal(DATA,model,params(k,:),paramsFixed)<limit
                     maxVal(maxVal<simVal)=simVal(maxVal<simVal);
                     minVal(minVal>simVal)=simVal(minVal>simVal);
-            if maxVal(2,:)>0.3
-                optParams = params(k,:);
-            end
-
+                end
             end
     end
         
diff --git a/scripts/getMaxMinPred.m b/scripts/getMaxMinPred.m
new file mode 100644
index 0000000000000000000000000000000000000000..e3feb055ac83a915bfb1085f4b30d5a77cd323a9
--- /dev/null
+++ b/scripts/getMaxMinPred.m
@@ -0,0 +1,36 @@
+function [boundries,bestParams,bestCost]=getMaxMinPred(simTime, model, params,paramsFixed,DATA)
+simVal=inf(2,length(simTime));
+
+minVal=simVal;
+maxVal=zeros(2,length(simTime));
+bestCost = 1e99;
+
+    for k=1:size(params,1)
+          try
+             [Gp,Ip,~,~,~,~] = simulateModelDalla(simTime, model, params(k,:),paramsFixed);
+             simVal=[Gp;Ip];
+          catch error
+              disp(error)
+          end
+
+            if k == 1
+                minVal=simVal;
+                maxVal=minVal;
+            end
+            if max(abs(simVal(:,:)))<inf  
+            costD = costFunDalla(model,params(k,:),paramsFixed);
+            if costD<79
+                    maxVal(maxVal<simVal)=simVal(maxVal<simVal);
+                    minVal(minVal>simVal)=simVal(minVal>simVal);   
+            cost = costFunDVal(model,params(k,:),paramsFixed,DATA);
+            if cost<bestCost
+                bestParams = params(k,:);
+                bestCost = cost;
+            end
+            end
+            end
+    end
+        
+boundries=table(maxVal,minVal,'VariableNames',{'Max','Min'},'RowNames',{'g','i'});
+
+end
\ No newline at end of file
diff --git a/scripts/get_maxmin.m b/scripts/get_maxmin.m
deleted file mode 100644
index bc4d77ce021a89998a36da46352b3a3c1db3d71c..0000000000000000000000000000000000000000
--- a/scripts/get_maxmin.m
+++ /dev/null
@@ -1,30 +0,0 @@
-function boundries=get_maxmin(simTime, model, params,paramsFixed)
-simVal=inf(2,length(simTime));
-
-minVal=simVal;
-maxVal=zeros(2,length(simTime));
-
-    for k=1:size(params,1)
-          try
-             [sim_muscle,sim_adipose,~,~] = simulateModel(simTime, model, params(k,:),paramsFixed);
-             simVal=[sim_muscle;sim_adipose];
-          catch error
-              disp(error)
-          end
-
-            if k == 1
-                minVal=simVal;
-                maxVal=minVal;
-            end
-            if max(abs(simVal(:,:)))<inf
-                %if abs(diff(simVal,1,2))<1e-4 
-                    maxVal(maxVal<simVal)=simVal(maxVal<simVal);
-                    minVal(minVal>simVal)=simVal(minVal>simVal);
-                %end
-
-            end
-    end
-        
-boundries=table(maxVal,minVal,'VariableNames',{'Max','Min'},'RowNames',{'m','a'});
-
-end
\ No newline at end of file
diff --git a/scripts/main.m b/scripts/main.m
index 91da64c0ec6b38fc3a0ab0328b8267bc150126e3..7fd5f7a6d26fc81bafbd0ebc8a114415761f919a 100644
--- a/scripts/main.m
+++ b/scripts/main.m
@@ -1,13 +1,14 @@
-% Do not change this script. 
 
 clear all
 rng('shuffle')
 
-Setup 
+% Setup 
 restoredefaultpath
+time = datestr(now,'yymmdd');
 
 %% Run optimizations
-parfor i = 1:100
+parfor i = 1:10
     seed=randi([0,32767],1);
-    EstimateParameters(seed,datestr(now,'yymmdd-HHMM')) % you can add more input if you want, but do not remove seed. 
+    EstimateParameters(seed,time)
+    EstimateParametersInputFat(seed,time)
 end
diff --git a/scripts/plotBrannmark.m b/scripts/plotBrannmark.m
deleted file mode 100644
index abe73583ec94f9a80c71445086abb08f36a1fec2..0000000000000000000000000000000000000000
--- a/scripts/plotBrannmark.m
+++ /dev/null
@@ -1,621 +0,0 @@
-%% Plot script
-restoredefaultpath
-addpath(genpath('./suport/SBPOP_PACKAGE_Rev_98'))
-addpath(genpath('./Models'))
-addpath(genpath('./Data'))
-
-wd = cd;
-try
-run('./suport/SBPOP_PACKAGE_Rev_98/installSBPOPpackageInitial.m')
-catch
-cd(wd);
-run('./suport/SBPOP_PACKAGE_Rev_98/installSBPOPpackage.m')
-end
-
-format shortg
-
-%% Load data
-wd = cd; % saves current directory
-cd data % Data directory
-load expdata_IR_dr.mat
-load expdata_IR_time.mat
-load expdata_IR_2p.mat
-load expdata_IR_int.mat
-load expdata_IRS1_dr.mat
-load expdata_IRS1_2p.mat
-load expdata_IRS1_time_3.mat
-load expdata_IRS1_time_30.mat
-load expdata_IRS1_time_180.mat
-load expdata_IRS1_ds.mat
-load expdata_IRS1307_dr.mat
-load expdata_IRS1307_time.mat
-load expdata_PKB308_dr.mat
-load expdata_PKB308_time.mat
-load expdata_PKB473_dr.mat
-load expdata_PKB473_time.mat
-load expdata_PKB473_2p.mat
-load expdata_PKB473_time_renorm.mat
-load expdata_AS160_dr.mat
-load expdata_AS160_time.mat
-load expdata_GLUCOSE_dr.mat
-load expdata_GLUCOSE_2p.mat
-load expdata_GLU_time.mat
-load expdata_S6K_time.mat
-load expdata_S6_time.mat
-load expdata_EC50.mat
-cd(wd) % returns to previous directory
-
-%% Load the model
-cd Models
-modelName = 'model_diabetes';
-optModel = SBmodel(strcat(modelName,'.txt'));
-tmpModel = optModel;
-SBPDmakeMEXmodel(optModel,modelName);
-cd ..
-
-%% Simulation options
-simOptions              = [];
-simOptions.method       = 'stiff';
-simOptions.maxnumsteps  = 1e6;
-simOptions.reltol       = 1e-12;
-simOptions.abstol       = 1e-9;
-simOptions.minstep      = 1e-10;
-
-[paramNames, param] = SBparameters(optModel);
-[stateNames,ODEs,initcond] = SBstates(optModel);
-[reactionNames] = SBreactions(optModel);
-figsize = [500 0 900 800];
-
-%% Defining parameter values
-
-size(param);
-tempParam = param';
-
-index_diabetes = ismember(paramNames,'diabetes')==1;
-index_ins = ismember(paramNames,'insulin')==1;
-index_gluc = ismember(paramNames,'gluc')==1;
-
-index_IR = ismember(stateNames,'IR')==1;
-index_GLUT4 = ismember(stateNames,'GLUT4')==1;
-index_cure_diabetes = ismember(paramNames,'cure_diabetes')==1;
-
-
-%% Baseline simulations
-% Parameters for Controls
-param_n = tempParam;                                % Parameter values
-param_n(index_diabetes)     = 1;                    % No diabetes feedback
-
-% Parameters for Diabetics
-param_d = tempParam;                                % Parameter values for diabetes
-param_d(index_diabetes) = 15.5/100;                 % Diabetes feedback 15.5 % of controls
-param_d(index_cure_diabetes) = 1;                   % Deacrese the dephosphorylation
-initcond_0_d = initcond;                            % Initial conditions
-initcond_0_d(index_IR) = 55;                        % Diabetes IR amount 55 % of controls
-initcond_0_d(index_GLUT4) = 50;                     % Diabetes GLUT4 amount 50 % of controls
-
-% Baseline simulations
-simdata_0_n = model_diabetes(0:50:500, initcond, param_n, simOptions);      % baseline simulation controls
-simdata_0_d = model_diabetes(0:50:500, initcond_0_d, param_d, simOptions);  % baseline simulation diabetes
- 
-%% Insulin stimulation simulations
-ins_conc = [0.01 0.03 0.1 0.3 1 10 100];
-ins_conc_highres = interp1([1:1:7],ins_conc,[1:0.1:7]);
-
-% Controls
-initcond_n = simdata_0_n.statevalues(end,:);    % using the steady state values from the baseline-simulation
-simdata_dr_n1a = model_diabetes(0:0.01:15, initcond_n, [param_n(1:end-1) ins_conc(1)], simOptions);   % 15 min simulation with 0.01 nM insulin
-simdata_dr_n1b = model_diabetes(0:0.01:15, initcond_n, [param_n(1:end-1) ins_conc(2)], simOptions);   % 15 min simulation with 0.03 nM insulin
-simdata_dr_n2 = model_diabetes(0:0.01:15, initcond_n, [param_n(1:end-1) ins_conc(3)], simOptions);    % 15 min simulation with 0.1 nM insulin
-simdata_dr_n3 = model_diabetes(0:0.01:15, initcond_n, [param_n(1:end-1) ins_conc(4)], simOptions);    % 15 min simulation with 0.3 nM insulin
-simdata_dr_n4 = model_diabetes(0:0.01:15, initcond_n, [param_n(1:end-1) ins_conc(5)], simOptions);    % 15 min simulation with 1 nM insulin
-simdata_60_n = model_diabetes(0:0.01:60, initcond_n, [param_n(1:end-1) ins_conc(6)], simOptions);     % 60 min simulation with 10 nM insulin
-simdata_180_n = model_diabetes(0:0.01:30, initcond_n, [param_n(1:end-1) ins_conc(7)], simOptions);   % 180 min simulation with 100 nM insulin
-simdata_ds1_n = model_diabetes(0:0.01:4, initcond_n, [param_n(1:end-1) 1.2], simOptions);       % 4 min simulation with 1.2 nM insulin for double step data
-initcond2_n = simdata_ds1_n.statevalues(end,:);                                 % using the statevalues from the ds1 simulation as initial conditions
-simdata_ds2_n = model_diabetes(0:0.01:4, initcond2_n, [param_n(1:end-1) 10], simOptions);        % 4 min simulation with 10 nM insulin for double step data
-
-% Diabetes
-initcond_d = simdata_0_d.statevalues(end,:);                                % using the steady state values from the baseline-simulation
-simdata_dr_d1a = model_diabetes(0:0.01:15, initcond_d, [param_d(1:end-1) ins_conc(1)], simOptions);   % 15 min simulation with 0.01 nM insulin   
-simdata_dr_d1b = model_diabetes(0:0.01:15, initcond_d, [param_d(1:end-1) ins_conc(2)], simOptions);   % 15 min simulation with 0.03 nM insulin
-simdata_dr_d2 = model_diabetes(0:0.01:15, initcond_d, [param_d(1:end-1) ins_conc(3)], simOptions);    % 15 min simulation with 0.1 nM insulin
-simdata_dr_d3 = model_diabetes(0:0.01:15, initcond_d, [param_d(1:end-1) ins_conc(4)], simOptions);    % 15 min simulation with 0.3 nM insulin
-simdata_dr_d4 = model_diabetes(0:0.01:15, initcond_d, [param_d(1:end-1) ins_conc(5)], simOptions);    % 15 min simulation with 1 nM insulin
-simdata_60_d = model_diabetes(0:0.01:60, initcond_d, [param_d(1:end-1) ins_conc(6)], simOptions);     % 60 min simulation with 10 nM insulin
-simdata_180_d = model_diabetes(0:0.01:30, initcond_d, [param_d(1:end-1) ins_conc(7)], simOptions);   % 180 min simulation with 100 nM insulin
-simdata_ds1_d = model_diabetes(0:0.01:4, initcond_d, [param_d(1:end-1) 1.2], simOptions);       % 4 min simulation with 1.2 nM insulin for double step data
-initcond2_d = simdata_ds1_d.statevalues(end,:);                                 % using the statevalues from the ds1 simulation as initial conditions
-simdata_ds2_d = model_diabetes(0:0.01:4, initcond2_d, [param_d(1:end-1) 10], simOptions);        % 4 min simulation with 10 nM insulin for double step data
-
-dose_n_highres = nan(length(ins_conc_highres)+1,13);
-dose_n_highres(1,:) = simdata_0_n.variablevalues(end,:);
-
-dose_d_highres = nan(length(ins_conc_highres)+1,13);
-dose_d_highres(1,:) = simdata_0_d.variablevalues(end,:);
-
-for k = 1:length(ins_conc_highres)
-simdata_n = model_diabetes(0:0.01:30, initcond_n, [param_n(1:end-1) ins_conc_highres(k)], simOptions);   % 180 min simulation with k nM insulin
-dose_n_highres(k+1,:) = simdata_n.variablevalues(1001,:);
-end
-
-for k = 1:length(ins_conc_highres)
-simdata_d = model_diabetes(0:0.01:30, initcond_d, [param_d(1:end-1) ins_conc_highres(k)], simOptions);   % 180 min simulation with k nM insulin
-dose_d_highres(k+1,:) = simdata_d.variablevalues(1001,:);
-end
-
-%% Glucose stimulation simulations
-% Controls
-param_n(index_gluc)=0.05; % Stimulating with 0.05 mM glucose
-simdata_dr_gluc_n0 = model_diabetes(0:0.01:30, initcond_n, [param_n(1:end-1) 0], simOptions);     % initial conditions from baseline simulation, 30 min simulation without insulin
-initcond_n1a = simdata_dr_n1a.statevalues(end,:);                                   % initial conditions from 0.01 nM insulin simulation (dr_n1a)
-simdata_dr_gluc_n1a = model_diabetes(0:0.01:30, initcond_n1a, [param_n(1:end-1) ins_conc(1)], simOptions);  % 30 min simulation with 0.01 nM insulin
-initcond_n1b = simdata_dr_n1b.statevalues(end,:);                                   % initial conditions from 0.03 nM insulin simulation (dr_n1b)
-simdata_dr_gluc_n1b = model_diabetes(0:0.01:30, initcond_n1b, [param_n(1:end-1) ins_conc(2)], simOptions);  % 30 min simulation with 0.03 nM insulin 
-initcond_n2 = simdata_dr_n2.statevalues(end,:);                                     % initial conditions from 0.1 nM insulin simulation (dr_n2)
-simdata_dr_gluc_n2 = model_diabetes(0:0.01:30, initcond_n2, [param_n(1:end-1) ins_conc(3)], simOptions);    % 30 min simulation with 0.1 nM insulin 
-initcond_n3 = simdata_dr_n3.statevalues(end,:);                                     % initial conditions from 0.3 nM insulin simulation (dr_n3)
-simdata_dr_gluc_n3 = model_diabetes(0:0.01:30, initcond_n3, [param_n(1:end-1) ins_conc(4)], simOptions);    % 30 min simulation with 0.3 nM insulin 
-initcond_n4 = simdata_dr_n4.statevalues(end,:);                                     % initial conditions from .01 nM insulin simulation (dr_n4)
-simdata_dr_gluc_n4 = model_diabetes(0:0.01:30, initcond_n4, [param_n(1:end-1) ins_conc(5)], simOptions);    % 30 min simulation with 1 nM insulin 
-initcond_n5 = simdata_60_n.statevalues(1501,:);                                      % initial conditions from 10 nM insulin simulation (60_n)
-simdata_dr_gluc_n5 = model_diabetes(0:0.01:30, initcond_n5, [param_n(1:end-1) ins_conc(6)], simOptions);    % 30 min simulation with 10 nM insulin 
-initcond_n6 = simdata_180_n.statevalues(1501,:);                                     % initial conditions from 100 nM insulin simulation (180_n)
-simdata_dr_gluc_n6 = model_diabetes(0:0.01:30, initcond_n6, [param_n(1:end-1) ins_conc(7)], simOptions);    % 30 min simulation with 100 nM insulin 
-
-% Diabetes
-param_d(index_gluc)=0.05; % Stimulating with 0.05 mM glucose
-simdata_dr_gluc_d0 = model_diabetes(0:0.01:30, initcond_d, [param_d(1:end-1) 0], simOptions);     % initial conditions from baseline simulation, 30 min simulation without insulin
-initcond_d1a = simdata_dr_d1a.statevalues(end,:);                                   % initial conditions from 0.01 nM insulin simulation (dr_d1a)
-simdata_dr_gluc_d1a = model_diabetes(0:0.01:30, initcond_d1a, [param_d(1:end-1) ins_conc(1)], simOptions);  % 30 min simulation with 0.01 nM insulin 
-initcond_d1b = simdata_dr_d1b.statevalues(end,:);                                   % initial conditions from 0.03 nM insulin simulation (dr_d1b)
-simdata_dr_gluc_d1b = model_diabetes(0:0.01:30, initcond_d1b, [param_d(1:end-1) ins_conc(2)], simOptions);  % 30 min simulation with 0.03 nM insulin 
-initcond_d2 = simdata_dr_d2.statevalues(end,:);                                     % initial conditions from 0.1 nM insulin simulation (dr_d2)
-simdata_dr_gluc_d2 = model_diabetes(0:0.01:30, initcond_d2, [param_d(1:end-1) ins_conc(3)], simOptions);    % 30 min simulation with 0.1 nM insulin 
-initcond_d3 = simdata_dr_d3.statevalues(end,:);                                     % initial conditions from 0.3 nM insulin simulation (dr_d3)
-simdata_dr_gluc_d3 = model_diabetes(0:0.01:30, initcond_d3, [param_d(1:end-1) ins_conc(4)], simOptions);    % 30 min simulation with 0.3 nM insulin 
-initcond_d4 = simdata_dr_d4.statevalues(end,:);                                     % initial conditions from .01 nM insulin simulation (dr_d4)
-simdata_dr_gluc_d4 = model_diabetes(0:0.01:30, initcond_d4, [param_d(1:end-1) ins_conc(5)], simOptions);    % 30 min simulation with 1 nM insulin 
-initcond_d5 = simdata_60_d.statevalues(1501,:);                                      % initial conditions from 10 nM insulin simulation (60_d)
-simdata_dr_gluc_d5 = model_diabetes(0:0.01:30, initcond_d5, [param_d(1:end-1) ins_conc(6)], simOptions);    % 30 min simulation with 10 nM insulin 
-initcond_d6 = simdata_180_d.statevalues(1501,:);                                     % initial conditions from 100 nM insulin simulation (180_d)
-simdata_dr_gluc_d6 = model_diabetes(0:0.01:30, initcond_d6, [param_d(1:end-1) ins_conc(7)], simOptions);    % 30 min simulation with 100 nM insulin 
-
-
-%% ExpData and SimData normalisations
-
-% Dose response
-% Collecting simulation data for variable values
-% Controls
-dose_n = [simdata_0_n.variablevalues(end,:);simdata_dr_n1a.variablevalues(1001,:); simdata_dr_n1b.variablevalues(1001,:); simdata_dr_n2.variablevalues(1001,:); simdata_dr_n3.variablevalues(1001,:); simdata_dr_n4.variablevalues(1001,:); simdata_60_n.variablevalues(1001,:); simdata_180_n.variablevalues(1001,:)];
-% Diabetes
-dose_d = [simdata_0_d.variablevalues(end,:);simdata_dr_d1a.variablevalues(1001,:); simdata_dr_d1b.variablevalues(1001,:); simdata_dr_d2.variablevalues(1001,:); simdata_dr_d3.variablevalues(1001,:); simdata_dr_d4.variablevalues(1001,:); simdata_60_d.variablevalues(1001,:); simdata_180_d.variablevalues(1001,:)];
-
-% Saving dose response data for IRp + IRip (variablevalue 1)
-%                               IRS1p (variablevalue 2)
-%                               IRS1307 (variablevalue 3)
-%                               PKB308 (variablevalue 4)
-%                               PKB473 (variablevalue 5)
-%                               AS160 (variablevalue 6)
-%                               S6 (variablevalue 9)
-% The doseresponse curves are normalized to 100 at maximum value
-% Baselines removed since data was treated that way
-% All dose response curves in this loop have insulin concentrations:
-% 1E-11,1E-10,3E-10,1E-9,1E-8,1E-7
-
-dr_n = zeros(length(ins_conc_highres),8);
-dr_d = zeros(length(ins_conc_highres),8);
-
-dr_variables = [1:6 9 7];
-
-%Controls
-for n = 1:8 
-dr_n(:,n) = [dose_n_highres(1,dr_variables(n)); dose_n_highres(3:end,dr_variables(n))];
-dr_n(:,n) = dr_n(:,n) - simdata_0_n.variablevalues(end,dr_variables(n)); %baseline removal
-dr_n(:,n) = dr_n(:,n)*100/max(dr_n(:,n));
-end
-
-%Diabetes
-for n = 1:8
-dr_d(:,n) = [dose_d_highres(1,dr_variables(n)); dose_d_highres(3:end,dr_variables(n))];
-dr_d(:,n) = dr_d(:,n) - simdata_0_d.variablevalues(end,dr_variables(n)); %baseline removal
-dr_d(:,n) = dr_d(:,n)*100/max(dr_d(:,n));
-end
-
-% Saving dose response data for GLUCOSE (variablevalue 11)
-% The doseresponse curve is normalized to 100 at its maximum value
-% Insulin concentrations: 1E-11,3E-11,1E-10,3E-10,1E-9,1E-8,1E-7
-% 
-% Collecting simulation data for glucose after insulin and glucose stimulation
-GLUCOSE_dr_n1 = [simdata_dr_gluc_n0.variablevalues(end,11) simdata_dr_gluc_n1a.variablevalues(end,11) simdata_dr_gluc_n1b.variablevalues(end,11) simdata_dr_gluc_n2.variablevalues(end,11) simdata_dr_gluc_n3.variablevalues(end,11) simdata_dr_gluc_n4.variablevalues(end,11) simdata_dr_gluc_n5.variablevalues(end,11) simdata_dr_gluc_n6.variablevalues(end,11)];
-GLUCOSE_dr_d1 = [simdata_dr_gluc_d0.variablevalues(end,11) simdata_dr_gluc_d1a.variablevalues(end,11) simdata_dr_gluc_d1b.variablevalues(end,11) simdata_dr_gluc_d2.variablevalues(end,11) simdata_dr_gluc_d3.variablevalues(end,11) simdata_dr_gluc_d4.variablevalues(end,11) simdata_dr_gluc_d5.variablevalues(end,11) simdata_dr_gluc_d6.variablevalues(end,11)];
-
-GLUCOSE_dr_n = GLUCOSE_dr_n1 - simdata_dr_gluc_n0.variablevalues(end,11);
-GLUCOSE_dr_n = GLUCOSE_dr_n*100/max(GLUCOSE_dr_n);
-GLUCOSE_dr_d = GLUCOSE_dr_d1 - simdata_dr_gluc_d0.variablevalues(end,11);
-GLUCOSE_dr_d = GLUCOSE_dr_d*100/max(GLUCOSE_dr_d);
-
-% Time curves
-
-% IRp + IRip
-% IR time has baseline subtracted and is normalised to max 100
-% 2p IR data is not normalised
-% Here, the data is combined to a 30 min dataset with baseline
-% Combination: Multiply expdata_time with expdata_2p without baseline divided by expdata_time at 10 min => expdata_time(1) = 0, expdata_time(9)=1. 
-%              Then add baseline from expdata_2p.
-comb_expdata_IR_time.response = expdata_IR_time.response*(expdata_IR_2p.response_n(2)-expdata_IR_2p.response_n(1))/expdata_IR_time.response(9); % normalising the 30 min data to the 10 min 2p data point
-comb_expdata_IR_time.response = comb_expdata_IR_time.response + expdata_IR_2p.response_n(1); % adding the baseline from the 2p data
-comb_expdata_IR_time.sem = expdata_IR_time.sem*(expdata_IR_2p.response_n(2)-expdata_IR_2p.response_n(1))/expdata_IR_time.response(9);
-IR_timep_n = simdata_180_n.variablevalues(1:3001,1);
-IR_time_n = IR_timep_n(round(expdata_IR_time.time*100)+1);
-IR_timep_d = simdata_180_d.variablevalues(1:3001,1);
-IR_time_d = IR_timep_d(round(expdata_IR_time.time*100)+1);
-scaleIR = lscov(IR_time_n,comb_expdata_IR_time.response); 
-
-% IRS1 
-% IRS1 time has baseline subtracted and is normalised to max 100
-% 2p IRS1 data is not normalised
-% Here, the data is combined to a 30 min dataset with baseline
-% Combination: Multiply expdata_time with expdata_2p without baseline divided by expdata_time at 10 min => expdata_time(1) = 0, expdata_time(9)=1. 
-%              Then add baseline from expdata_2p.
-comb_expdata_IRS1_time.response = expdata_IRS1_time_30.response*(expdata_IRS1_2p.response_n(2)-expdata_IRS1_2p.response_n(1))/expdata_IRS1_time_30.response(9); % normalising the 30 min data to the 10 min 2p data point
-comb_expdata_IRS1_time.response = comb_expdata_IRS1_time.response + expdata_IRS1_2p.response_n(1); % adding the baseline from the 2p data
-comb_expdata_IRS1_time.sem = expdata_IRS1_time_30.sem*(expdata_IRS1_2p.response_n(2)-expdata_IRS1_2p.response_n(1))/expdata_IRS1_time_30.response(9);
-IRS1_timep_n = simdata_180_n.variablevalues(1:3001,2);
-IRS1_time_n = IRS1_timep_n(round(expdata_IRS1_time_30.time*100)+1);
-IRS1_timep_d = simdata_180_d.variablevalues(1:3001,2);
-IRS1_time_d = IRS1_timep_d(round(expdata_IRS1_time_30.time*100)+1);
-scaleIRS1 = lscov(IRS1_time_n,comb_expdata_IRS1_time.response); 
-
-% IRS1p 10nm insulin 3 min 
-% baseline not subtracted
-% data normalised to max 100, simdata scaled with lscov
-IRS1_time3p = simdata_60_n.variablevalues(1:301,2);
-IRS1_time3 = IRS1_time3p(round(expdata_IRS1_time_3.time.*100)+1);
-scaleIRS1_10 = lscov(IRS1_time3,expdata_IRS1_time_3.response);
-
-% IRS1p double step (1.2 + 10 nm)
-% baseline not subtracted
-% simdata scaled with lscov
-IRS1_dsp = [simdata_ds1_n.variablevalues(:,2); simdata_ds2_n.variablevalues(:,2)];
-IRS1_ds = IRS1_dsp(round(expdata_IRS1_ds.time.*100)+1); 
-scaleIRS1ds = lscov(IRS1_ds,expdata_IRS1_ds.response);
-
-% IRS1307
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, simulation scaled with lscov
-IRS1307_time_n = simdata_60_n.variablevalues(:,3);
-IRS1307_time_n = IRS1307_time_n(round(expdata_IRS1307_time.time*100)+1); %choosing correct values
-IRS1307_time_d = simdata_60_d.variablevalues(:,3);
-IRS1307_time_d = IRS1307_time_d(round(expdata_IRS1307_time.time*100)+1);
-scaleIRS1307 = lscov(IRS1307_time_n,expdata_IRS1307_time.response_n); %scale parameter for IRS1307
-
-% PKB308 
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, simulation scaled with lscov
-PKB308_time_n = simdata_60_n.variablevalues(:,4);
-PKB308_time_n = PKB308_time_n(round(expdata_PKB308_time.time*100+1));
-PKB308_time_d = simdata_60_d.variablevalues(:,4);
-PKB308_time_d = PKB308_time_d(round(expdata_PKB308_time.time*100+1));
-scalePKB308 = lscov(PKB308_time_n,expdata_PKB308_time.response_n);
-
-% PKB473
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, ExpData normalised to max 100. Normalised to 10 point in 2p ExpData to get relationship between n and d
-% simulation scaled with lscov
-% Combination: Multiply expdata_time with expdata_2p divided by expdata_time at 10 min => expdata_time(9)=1. 
-comb_expdata_PKB473_time.response_n = expdata_PKB473_time_renorm.response_n*expdata_PKB473_2p.response_n(2); % normalising the 60 min data to the 10 min 2p data point
-comb_expdata_PKB473_time.sem_n = expdata_PKB473_time_renorm.sem_n*expdata_PKB473_2p.response_n(2);
-comb_expdata_PKB473_time.response_d = expdata_PKB473_time_renorm.response_d*expdata_PKB473_2p.response_d(2); % normalising the 60 min data to the 10 min 2p data point
-comb_expdata_PKB473_time.sem_d = expdata_PKB473_time_renorm.sem_d*expdata_PKB473_2p.response_d(2);
-
-PKB473_time_n = simdata_60_n.variablevalues(:,5);
-PKB473_time_n = PKB473_time_n(round(expdata_PKB473_time_renorm.time*100+1));
-PKB473_time_d = simdata_60_d.variablevalues(:,5);
-PKB473_time_d = PKB473_time_d(round(expdata_PKB473_time_renorm.time*100+1));
-scalePKB473 = lscov(PKB473_time_n,comb_expdata_PKB473_time.response_n);
-
-% AS160
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, simulation scaled with lscov
-AS160_time_n = simdata_60_n.variablevalues(:,6);
-AS160_time_n = AS160_time_n(round(expdata_AS160_time.time*100+1));
-AS160_time_d = simdata_60_d.variablevalues(:,6);
-AS160_time_d = AS160_time_d(round(expdata_AS160_time.time*100+1));
-scaleAS160 = lscov(AS160_time_n,expdata_AS160_time.response_n);
-
-% Glucose uptake
-% 0 and 100 nM insulin 15 min followed by 0.05 nM glucose 30 min
-GLUCOSE_n = [simdata_dr_gluc_n0.variablevalues(end,11); simdata_dr_gluc_n6.variablevalues(end,11)];
-GLUCOSE_d = [simdata_dr_gluc_d0.variablevalues(end,11); simdata_dr_gluc_d6.variablevalues(end,11)];
-scaleGLUCOSE = lscov(GLUCOSE_n,expdata_GLUCOSE_2p.response_n);
-
-% S6K 
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, simulation scaled with lscov
-S6K_time_n = simdata_60_n.variablevalues(:,8);
-S6K_time_n = S6K_time_n(round(expdata_S6K_time.time*100)+1);
-S6K_time_d = simdata_60_d.variablevalues(:,8);
-S6K_time_d = S6K_time_d(round(expdata_S6K_time.time*100)+1);
-scaleS6K = lscov(S6K_time_n, expdata_S6K_time.response_n);
-
-% S6
-% 10nm insulin 60min normal (n) and diabetes (d)
-% baseline not subtracted, simulation scaled with lscov
-S6_time_n = simdata_60_n.variablevalues(:,9);
-S6_time_n = S6_time_n(round(expdata_S6_time.time*100)+1);
-S6_time_d = simdata_60_d.variablevalues(:,9);
-S6_time_d = S6_time_d(round(expdata_S6_time.time*100)+1);
-scaleS6 = lscov(S6_time_n, expdata_S6_time.response_n);
-
-%% Plot
-time10 = 0:0.01:10;
-time15 = 0:0.01:15;
-time30 = 0:0.01:30;
-time60 = 0:0.01:60;
-conc = 1e-9*[1e-3 ins_conc([1 3:end])];
-conc_highres = interp1([1:1:7],conc,[1:0.1:7]);
-
-opengl hardware
-set(0,'DefaultLineLineSmoothing','on');
-set(0,'DefaultPatchLineSmoothing','on');
-opengl('OpenGLLineSmoothingBug',1);
-
-set(figure(11),'Position', figsize);
-figure(11)
-% Dose response
-subplot(5,4,1)
-semilogx(conc_highres, dr_n(:,1)','b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,1)','r-', 'Linewidth',2)
-errorbar(conc,expdata_IR_dr.response_n,expdata_IR_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc,expdata_IR_dr.response_d,expdata_IR_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8)
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('IR-YP, % of max','FontSize',8)
-text(2e-12,95,'IR-YP','FontSize',9)
-
-subplot(5,4,2)
-errorbar(expdata_IR_time.time,comb_expdata_IR_time.response,comb_expdata_IR_time.sem,'bo','MarkerFaceColor','b','MarkerSize',6)
-hold on
-errorbar([0 10],expdata_IR_2p.response_d,expdata_IR_2p.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-plot(time30,simdata_180_n.variablevalues(1:3001,1)*scaleIR,'b-', 'Linewidth',2)
-plot(time30,simdata_180_d.variablevalues(1:3001,1)*scaleIR,'r-', 'Linewidth',2)
-axis([0 30 0 2.5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8)
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('IR-YP, a.u.','FontSize',8)
-text(15,0.95*2.5,'IR-YP','FontSize',9)
-
-subplot(5,4,3)
-errorbar(expdata_IR_int.time,expdata_IR_int.response, expdata_IR_int.sem,'bo','MarkerFaceColor','b','MarkerSize',6)
-hold on
-plot(time30,simdata_60_n.variablevalues(1:3001,12)*100,'b-', 'Linewidth',2)
-plot(time30,simdata_60_d.variablevalues(1:3001,12)*100,'r-', 'Linewidth',2)
-axis([0 30 0 4])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8)
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('Internalized IR, a.u.','FontSize',8)
-text(1.5,0.95*4,'Internalized IR','FontSize',9)
-
-subplot(5,4,4)
-errorbar(expdata_IRS1_ds.time,expdata_IRS1_ds.response,expdata_IRS1_ds.sem,'bo','MarkerFaceColor','b','MarkerSize',6)
-hold on
-plot([simdata_ds1_n.time 4+simdata_ds2_n.time]',scaleIRS1ds*[simdata_ds1_n.variablevalues(:,2); simdata_ds2_n.variablevalues(:,2)],'b-', 'Linewidth',2)
-plot([simdata_ds1_d.time 4+simdata_ds2_d.time]',scaleIRS1ds*[simdata_ds1_d.variablevalues(:,2); simdata_ds2_d.variablevalues(:,2)],'r-', 'Linewidth',2)
-axis([0 8 0 125])
-set(gca,'YTick',[0 50 100],'XTick',[0 4 8],'FontSize',8)
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('IRS1-YP, a.u.','FontSize',8)
-text(4,0.95*125,'IRS1-YP','FontSize',9)
-
-subplot(5,4,5)
-semilogx(conc_highres, dr_n(:,2)','b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,2)','r-', 'Linewidth',2)
-errorbar(conc,expdata_IRS1_dr.response_n,expdata_IRS1_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc,expdata_IRS1_dr.response_d,expdata_IRS1_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8)
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('IRS1-YP, % of max','FontSize',8)
-text(2e-12,95,'IRS1-YP','FontSize',9)
-
-subplot(5,4,6)
-errorbar(expdata_IRS1_time_30.time,comb_expdata_IRS1_time.response,comb_expdata_IRS1_time.sem,'bo','MarkerFaceColor','b','MarkerSize',6)
-hold on
-errorbar([0 10],expdata_IRS1_2p.response_d,expdata_IRS1_2p.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-plot(time30,scaleIRS1*simdata_180_n.variablevalues(1:3001,2),'b-', 'Linewidth',2)
-plot(time30,scaleIRS1*simdata_180_d.variablevalues(1:3001,2),'r-', 'Linewidth',2)
-axis([0 30 0 5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8)
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('IRS1-YP, a.u.','FontSize',8)
-text(15,0.95*5,'IRS1-YP','FontSize',9)
-
-subplot(5,4,7)
-semilogx(conc_highres, dr_n(:,3)','b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,3)','r-', 'Linewidth',2)
-errorbar(conc,expdata_IRS1307_dr.response_n,expdata_IRS1307_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc,expdata_IRS1307_dr.response_d,expdata_IRS1307_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8)
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('IRS1-S307P, % of max','FontSize',8)
-text(2e-12,95,'IRS1-S307P','FontSize',9)
-
-subplot(5,4,8)
-plot(time30,scaleIRS1307*simdata_60_n.variablevalues(1:3001,3),'b-', 'Linewidth',2)
-hold on
-plot(time30,scaleIRS1307*simdata_60_d.variablevalues(1:3001,3),'r-', 'Linewidth',2)
-errorbar(expdata_IRS1307_time.time,expdata_IRS1307_time.response_n,expdata_IRS1307_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_IRS1307_time.time,expdata_IRS1307_time.response_d,expdata_IRS1307_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([0 30 0 5.5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('IRS1-S307P, % of max','FontSize',8)
-text(15,0.95*5.5,'IRS1-S307P','FontSize',9)
-
-subplot(5,4,9)
-semilogx(conc_highres, dr_n(:,4)', 'b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,4)', 'r-', 'Linewidth',2)
-errorbar(conc, expdata_PKB308_dr.response_n,expdata_PKB308_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc, expdata_PKB308_dr.response_d,expdata_PKB308_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8);
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('PKB-T308P, % of max','FontSize',8)
-text(2e-12,95,'PKB-T308P','FontSize',9)
-
-subplot(5,4,10)
-plot(time30,scalePKB308.*simdata_60_n.variablevalues(1:3001,4),'b-', 'Linewidth',2)
-hold on
-plot(time30,scalePKB308.*simdata_60_d.variablevalues(1:3001,4),'r-', 'Linewidth',2)
-errorbar(expdata_PKB308_time.time,expdata_PKB308_time.response_n,expdata_PKB308_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_PKB308_time.time,expdata_PKB308_time.response_d,expdata_PKB308_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([0 30 0 3.5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('PKB-T308P, a.u.','FontSize',8)
-text(15,0.95*3.5,'PKB-T308P','FontSize',9)
-
-subplot(5,4,11)
-semilogx(conc_highres, dr_n(:,5)', 'b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,5)', 'r-', 'Linewidth',2)
-errorbar(conc, expdata_PKB473_dr.response_n,expdata_PKB473_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc, expdata_PKB473_dr.response_d,expdata_PKB473_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8);
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('PKB-S473P, % of max','FontSize',8)
-text(2e-12,95,'PKB-S473P','FontSize',9)
-
-subplot(5,4,12)
-plot(time30,scalePKB473*simdata_60_n.variablevalues(1:3001,5),'b-', 'Linewidth',2)
-hold on
-plot(time30,scalePKB473*simdata_60_d.variablevalues(1:3001,5),'r-', 'Linewidth',2)
-errorbar(expdata_PKB473_time_renorm.time,comb_expdata_PKB473_time.response_n,comb_expdata_PKB473_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_PKB473_time_renorm.time,comb_expdata_PKB473_time.response_d,comb_expdata_PKB473_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([0 30 0 2.5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('PKB-S473P, a.u.','FontSize',8)
-text(15,0.95*2.5,'PKB-S473P','FontSize',9)
-
-subplot(5,4,15)
-plot(time30,simdata_60_n.variablevalues(1:3001,7),'b-', 'Linewidth',2)
-hold on
-plot(time30,simdata_60_d.variablevalues(1:3001,7),'r-', 'Linewidth',2)
-axis([0 30 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('Active mTORC1, a.u.','FontSize',8)
-text(2,0.95*100,'Active mTORC1','FontSize',9)
-
-subplot(5,4,13)
-semilogx(conc_highres, dr_n(:,6)', 'b-', 'Linewidth',2)
-hold on
-plot(conc_highres, dr_d(:,6)', 'r-', 'Linewidth',2)
-errorbar(conc, expdata_AS160_dr.response_n,expdata_AS160_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(conc, expdata_AS160_dr.response_d,expdata_AS160_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8);
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('AS160-T642P, % of max','FontSize',8)
-text(2e-12,90,{'AS160-';'T642P'},'FontSize',9)
-
-subplot(5,4,14)
-plot(time30,scaleAS160*simdata_60_n.variablevalues(1:3001,6),'b-', 'Linewidth',2)
-hold on
-plot(time30,scaleAS160*simdata_60_d.variablevalues(1:3001,6),'r-', 'Linewidth',2)
-errorbar(expdata_AS160_time.time,expdata_AS160_time.response_n,expdata_AS160_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_AS160_time.time,expdata_AS160_time.response_d,expdata_AS160_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([0 30 0 3])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('AS160-T642P, a.u.','FontSize',8)
-text(15,0.95*3,'AS160-T642P','FontSize',9)
-
-subplot(5,4,17)
-semilogx(expdata_GLUCOSE_dr.conc,GLUCOSE_dr_n,'b-', 'Linewidth',2)
-hold on
-plot(expdata_GLUCOSE_dr.conc,GLUCOSE_dr_d,'r-', 'Linewidth',2)
-errorbar(expdata_GLUCOSE_dr.conc,expdata_GLUCOSE_dr.response_n,expdata_GLUCOSE_dr.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_GLUCOSE_dr.conc,expdata_GLUCOSE_dr.response_d,expdata_GLUCOSE_dr.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6)
-axis([1e-12 1e-7 0 100])
-set(gca,'YTick',[0 50 100],'XTick',[1e-11 1e-10 1e-9 1e-8 1e-7],'FontSize',8);
-box off
-xlabel('[Insulin], M','FontSize',8)
-ylabel('Glucose uptake, % of max','FontSize',8)
-text(2e-12,90,{'Glucose';'uptake'},'FontSize',9)
-
-subplot(5,4,18)
-errorbar(1,1e-3*expdata_GLUCOSE_2p.response_n(1),1e-3*expdata_GLUCOSE_2p.sem_n(1),'k')
-hold on
-errorbar(3,1e-3*expdata_GLUCOSE_2p.response_n(2),1e-3*expdata_GLUCOSE_2p.sem_n(2),'k')
-bar([1 3 2 4],[1e-3*expdata_GLUCOSE_2p.response_n;1e-3*scaleGLUCOSE*GLUCOSE_n],'b')
-errorbar(6,1e-3*expdata_GLUCOSE_2p.response_d(1),1e-3*expdata_GLUCOSE_2p.sem_d(2),'k')
-errorbar(8,1e-3*expdata_GLUCOSE_2p.response_d(2),1e-3*expdata_GLUCOSE_2p.sem_d(2),'k')
-bar([6 8 7 9],[1e-3*expdata_GLUCOSE_2p.response_d;1e-3*scaleGLUCOSE*GLUCOSE_d],'r')
-xlabel({'-       +           -       +';'Insulin'},'FontSize',8)
-box off
-ylabel('Glucose uptake, cpm*1e-3','FontSize',8)
-set(gca,'YTick',[0 1 2 3 4],'XTick',[],'FontSize',8);
-text(1,0.95*4,'Glucose uptake','FontSize',9)
-
-subplot(5,4,16)
-plot(time30,scaleS6K*simdata_60_n.variablevalues(1:3001,8),'b-', 'Linewidth',2)
-hold on
-plot(time30,scaleS6K*simdata_60_d.variablevalues(1:3001,8),'r-', 'Linewidth',2)   
-errorbar(expdata_S6K_time.time,expdata_S6K_time.response_n,expdata_S6K_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)        
-errorbar(expdata_S6K_time.time,expdata_S6K_time.response_d,expdata_S6K_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6')  
-axis([0 30 0 3.5])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('S6K-T389P, a.u.','FontSize',8)
-text(15,0.95*3.5,'S6K-T389P','FontSize',9)
-
-subplot(5,4,19)
-plot(time30,simdata_60_n.variablevalues(1:3001,13),'b-', 'Linewidth',2)
-hold on
-plot(time30,simdata_60_d.variablevalues(1:3001,13),'r-', 'Linewidth',2)
-axis([0 30 0 2.5])
-set(gca,'YTick',[0 1 2],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('Active mTORC2, a.u.','FontSize',8)
-text(5,0.95*2.5,'Active mTORC2','FontSize',9)
-
-subplot(5,4,20)
-plot(time30,scaleS6*simdata_60_n.variablevalues(1:3001,9),'b-', 'Linewidth',2) 
-hold on
-plot(time30,scaleS6*simdata_60_d.variablevalues(1:3001,9),'r-', 'Linewidth',2) 
-errorbar(expdata_S6_time.time,expdata_S6_time.response_n,expdata_S6_time.sem_n,'bo','MarkerFaceColor','b','MarkerSize',6)
-errorbar(expdata_S6_time.time,expdata_S6_time.response_d,expdata_S6_time.sem_d,'ro','MarkerFaceColor','r','MarkerSize',6) 
-axis([0 30 0 4])
-set(gca,'YTick',[0 1 2 3 4 5],'XTick',[0 10 20 30],'FontSize',8);
-box off
-xlabel('Time, min','FontSize',8)
-ylabel('S6-S235/236P, a.u.','FontSize',8)
-text(15,0.9*4,{'S6-S235/','236P'},'FontSize',9)
-
-
diff --git a/scripts/plotDalla.m b/scripts/plotDalla.m
deleted file mode 100644
index 61e886df5d1aa1550659fb369b5bf1e6690abe67..0000000000000000000000000000000000000000
--- a/scripts/plotDalla.m
+++ /dev/null
@@ -1,110 +0,0 @@
-% close all
-function []=plotDalla(sim_time,model_name,optParam)
-
-% load data and bounderies
-
-load('EGP_ori.mat');
-load('EGP_dalla.mat');
-load('EGP_high.mat');
-load('EGP_low.mat');
-
-load('S_ori.mat');
-load('S_dalla.mat');
-load('S_high.mat');
-load('S_low.mat');
-
-load('pi_ori.mat');
-load('pi_dalla.mat');
-load('pi_high.mat');
-load('pi_low.mat');
-
-load('pg_ori.mat');
-load('pg_dalla.mat');
-load('pg_high.mat');
-load('pg_low.mat');
-
-load('Ra_ori.mat');
-load('Ra_dalla.mat');
-load('Ra_high.mat');
-load('Ra_low.mat');
-
-load('GU_ori.mat');
-load('GU_dalla.mat');
-load('GU_high.mat');
-load('GU_low.mat');
-
-% simulate
-sim_time = [pg_dalla(1,1):0.01:pg_dalla(end,1)];
-paramsAll=[paramsFixed optParam];
-[inits] = simulateSteadyState(model, paramsAll);
-simModel = model(sim_time, inits, paramsAll); 
-
-% plot
-figure()
-set(0,'DefaultLineLineWidth',3)
-
-subplot(3,2,1)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'G')),'k-');
-hold on
-plot(pg_dalla(:,1),pg_dalla(:,2),'k.')
-hold on
-plotUncertaintyD(pg_high,pg_low,sim_time)
-ylabel('(mg/dl)')
-title('Plasma Glucose')
-
-subplot(3,2,2)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'I')),'k-');
-hold on
-plot(pi_dalla(:,1),pi_dalla(:,2),'k.')
-hold on
-plotUncertaintyD(pi_high,pi_low,sim_time)
-ylabel('(pmol/L)')
-title('Plasma Insulin')
-
-subplot(3,2,3)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'EGP')),'k-');
-hold on
-plot(EGP_dalla(:,1),EGP_dalla(:,2),'k.')
-hold on 
-plotUncertaintyD(EGP_high,EGP_low,sim_time)
-hold off
-ylabel('(mg/kg/min)')
-title('EGP')
-
-subplot(3,2,4)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'Ra')),'k-');
-hold on
-plot(Ra_dalla(:,1),Ra_dalla(:,2),'k.')
-hold on
-plotUncertaintyD(Ra_high,Ra_low,sim_time)
-ylabel('(mg/kg/min)')
-title('Glucose Rate of Appearance')
-
-subplot(3,2,5)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'U')),'k-');
-hold on
-plot(GU_dalla(:,1),GU_dalla(:,2),'k.')
-hold on
-plotUncertaintyD(GU_high,GU_low,sim_time)
-title('Glucose utilization')
-
-subplot(3,2,6)
-plot(simModel.time, simModel.variablevalues(:,ismember(IQMvariables(model),'S')),'k-');
-hold on
-plot(S_dalla(:,1),S_dalla(:,2),'k.')
-hold on
-plotUncertaintyD(S_high,S_low,sim_time)
-ylabel('(pmol/kg/min)')
-title('Insulin secretion S')
-
-
-function [] = plotUncertaintyD(high,low,sim_time)
-    high=interp1q(high(:,1),high(:,2),sim_time')';
-    low=interp1q(low(:,1),low(:,2),sim_time')';
-    h=fill([sim_time fliplr(sim_time)], [high fliplr(low)],'k','FaceAlpha',.3,'EdgeAlpha',.3);
-    set(h,'EdgeColor','None');
-    xlabel('time (min)');
-    ylabel('mg/kg/min');
-end
-
-end
\ No newline at end of file
diff --git a/scripts/results/GlucoseInsulin.png b/scripts/results/GlucoseInsulin.png
new file mode 100644
index 0000000000000000000000000000000000000000..91c37f98f3fa49396ada3f46f016a5d10cd528e1
Binary files /dev/null and b/scripts/results/GlucoseInsulin.png differ
diff --git a/scripts/results/M4m_allGoodParams.mat b/scripts/results/M4m_allGoodParams.mat
deleted file mode 100644
index ef6618720bfd55d11a34c117957a4e283a658013..0000000000000000000000000000000000000000
Binary files a/scripts/results/M4m_allGoodParams.mat and /dev/null differ
diff --git a/scripts/results/M4m_opt(28.86930).mat b/scripts/results/M4m_opt(28.86930).mat
deleted file mode 100644
index c8bcdbd0ba9c4bb72db7170c4441c040b10e8e81..0000000000000000000000000000000000000000
Binary files a/scripts/results/M4m_opt(28.86930).mat and /dev/null differ
diff --git a/scripts/results/allAcceptableParams.mat b/scripts/results/allAcceptableParams.mat
index 7d55b9dde43e2e09e890192f74c57c86634dfccb..752df77481d63422006968d3cca11626f1dec551 100644
Binary files a/scripts/results/allAcceptableParams.mat and b/scripts/results/allAcceptableParams.mat differ
diff --git a/scripts/results/boundriesVal.mat b/scripts/results/boundriesVal.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1278fb5064e5923eeb7f155814e6d95b0f38bd55
Binary files /dev/null and b/scripts/results/boundriesVal.mat differ
diff --git a/scripts/results/opt(12.92335).mat b/scripts/results/opt(12.92335).mat
deleted file mode 100644
index bba2ba9a18534124f99e84b35807718ac52fcbc2..0000000000000000000000000000000000000000
Binary files a/scripts/results/opt(12.92335).mat and /dev/null differ
diff --git a/scripts/results/optParam(24.549).mat b/scripts/results/optParam(24.549).mat
new file mode 100644
index 0000000000000000000000000000000000000000..21efda022b28a20f0ef381857899825b9d64374d
Binary files /dev/null and b/scripts/results/optParam(24.549).mat differ
diff --git a/scripts/results/optParamF.mat b/scripts/results/optParamF.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5ab8b57821c6f9d8a3c0b1b0a9f2086b17d9cfd2
Binary files /dev/null and b/scripts/results/optParamF.mat differ
diff --git a/scripts/results/optParamsPred.mat b/scripts/results/optParamsPred.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a12354a1141fe2799d3bccec197328e44dfa432f
Binary files /dev/null and b/scripts/results/optParamsPred.mat differ
diff --git a/scripts/runParallel.sh b/scripts/runParallel.sh
index 0e9b29861b8fc0f7d4b2b88d35d7e8b4619f0be6..76c6b8335e3adc8e8c627905098e126cf817afda 100644
--- a/scripts/runParallel.sh
+++ b/scripts/runParallel.sh
@@ -1,22 +1,27 @@
 #!/bin/bash
 #
-# Use 2 nodes, i.e. 2x32 cores, for 1 hour
+# Use 2 nodes, i.e. 2x32? cores, for 1 hour
 #SBATCH -N 2 
 #SBATCH -t 02:00:00
 
+# Update the project if possible. 
+#git pull
+
 # Load the Matlab module
 module load MATLAB/R2018b-nsc1 
 module add gcc/8.2.0 
 
 MATLAB='matlab -nodesktop -nodisplay -singleCompThread'
-
-job=EstimateParameters
+# The name of the Matlab script (without .m)
+job=EstimateParametersInputFat
 now=$(date '+%Y%m%d-%H%M%S')
 
+# Setup things first. 
 ${MATLAB} -r "Setup(); exit"
 
+# Run main task. Note the explicit "exit" to exit Matlab.
 for i in $(seq 1 256);do
-  srun -N1 -n1  ${MATLAB} -r "EstimateParametersInput($RANDOM,'${now}');exit" > ./Log/${job}_${i}.log &
+  srun -N1 -n1  ${MATLAB} -r "${job}($RANDOM,'${now}');exit" > ./Log/${job}_${i}.log &
 done
 wait # needed to prevent the script from exiting before all Matlab tasks are done
 
diff --git a/scripts/sampleParams.m b/scripts/sampleParams.m
index 65dfaa69d61991f1bdd0e0982b65d8a59a0ab876..09f446eb541d7518db89a8f330a6ef7d44bf6122 100644
--- a/scripts/sampleParams.m
+++ b/scripts/sampleParams.m
@@ -4,10 +4,12 @@ function [sampledParams] = sampleParams(nsets,allParams)
 % allocation
 np=size(allParams,2)-1;
 sets=floor(linspace(1,size(allParams,1),nsets));
-sampledParams=nan(nsets,np);
+sampledParams=[]; %nan(nsets*np,np);
 
 % sample
-allParams=sortrows(allParams,(np+1));
-sampledParams(:,:)=allParams(sets,1:np);
+for p = 1:np+1
+allParams = sortrows(allParams,p);
+sampledParams = [sampledParams;allParams(sets,1:np)];
+end
 
 end
diff --git a/scripts/sample_params.m b/scripts/sample_params.m
deleted file mode 100644
index b91e87f51e8f8aeb206e0ccc60c24d6e825c74d5..0000000000000000000000000000000000000000
--- a/scripts/sample_params.m
+++ /dev/null
@@ -1,13 +0,0 @@
-
-function [sampledParams] = sample_params(nsets,allParams)
-
-% allocation
-np=size(allParams,2)-1;
-sets=floor(linspace(1,size(allParams,1),nsets));
-sampledParams=nan(nsets,np);
-
-% sample
-allParams=sortrows(allParams,(np+1));
-sampledParams(:,:)=allParams(sets,1:np);
-
-end
diff --git a/scripts/simInput.m b/scripts/simInput.m
new file mode 100644
index 0000000000000000000000000000000000000000..ffecbab8543ec5355928ee32100f1d7c1590de5f
--- /dev/null
+++ b/scripts/simInput.m
@@ -0,0 +1,56 @@
+
+% This script was used to simulate the input curves of EGP, Ra, INS, GLUT4m
+% for simulation of sub module M4m.
+
+load ('Coppack.mat'); 
+% 
+model_name='M0';
+model=str2func(model_name);
+
+[pNames,paramsMod]=IQMparameters(model); 
+[inits] = simulateSteadyStated(model, paramsMod);
+sim = model(DATA.time, inits, paramsMod');  
+
+EGP = sim.variablevalues(:,ismember(IQMvariables(model),'EGP'))'; 
+Ra = sim.variablevalues(:,ismember(IQMvariables(model),'Ra'))';
+INS = sim.statevalues(:,ismember(IQMstates(model),'INS'))'; 
+GLUT4m = sim.statevalues(:,ismember(IQMstates(model),'GLUT4m'))';
+G_t = sim.statevalues(:,ismember(IQMstates(model),'G_t'))';  
+G_p = sim.statevalues(:,ismember(IQMstates(model),'G_p'))'; 
+I = sim.variablevalues(:,ismember(IQMvariables(model),'I'))'; 
+
+figure()
+subplot(2,2,1); hold on
+plot(sim.time, EGP)
+title('EGP')
+subplot(2,2,2); hold on
+plot(sim.time, Ra)
+title('Ra')
+subplot(2,2,3); hold on
+plot(sim.time, INS)
+title('INS')
+subplot(2,2,4); hold on
+plot(sim.time, GLUT4m)
+title('GLUT4m')
+
+figure()
+subplot(2,2,1); hold on
+plot(sim.time, G_t)
+title('G_t')
+subplot(2,2,2); hold on
+plot(sim.time, G_p)
+title('G_p')
+subplot(2,2,3); hold on
+plot(sim.time, I)
+title('I')
+
+
+
+save('./data/EGP.mat','EGP')
+save('./data/Ra.mat','Ra')
+save('./data/INS.mat','INS')
+save('./data/GLUT4m.mat','GLUT4m')
+save('./data/G_p.mat','G_p')
+save('./data/G_t.mat','G_t')
+save('./data/I.mat','I')
+
diff --git a/scripts/simulateModel.m b/scripts/simulateModel.m
index ae050ad582649be2b7b3d262fdc8cb2464e310a2..ef64d92600bec6923cdf1b136cd8b21ccb98ccdb 100644
--- a/scripts/simulateModel.m
+++ b/scripts/simulateModel.m
@@ -2,12 +2,12 @@
 function [sim_muscle,sim_adipose,sim_liver,sim_ii] = simulateModel(time, model, params,paramsFixed)
  
    paramsAll=[paramsFixed params];
-   [inits] = simulateSteadyState(model, paramsAll);
-   sim = model(time, inits, paramsAll'); 
+   [inits] = simulateSteadyStated(model, paramsAll);
+   sim = model(time, inits, paramsAll); 
    sim_muscle = sim.variablevalues(:,ismember(IQMvariables(model),'U_idm'))'; 
    sim_adipose = sim.variablevalues(:,ismember(IQMvariables(model),'U_idf'))'; 
    sim_liver = sim.variablevalues(:,ismember(IQMvariables(model),'U_idl'))';
    
-   sim_ii = paramsAll(end).*ones(size(sim_muscle));
+   sim_ii = params(1).*ones(size(sim_muscle));
 
 end
diff --git a/scripts/simulateModelDalla.m b/scripts/simulateModelDalla.m
new file mode 100644
index 0000000000000000000000000000000000000000..8589635947d1cf35cce3d7f2280af11711b93e4e
--- /dev/null
+++ b/scripts/simulateModelDalla.m
@@ -0,0 +1,14 @@
+
+function [Gp,Ip,EGP,Ra,U_id,S] = simulateModelDalla(time, model, params,paramsFixed)
+
+   paramsAll=[paramsFixed params];
+   [inits] = simulateSteadyStated(model, paramsAll);
+   sim = model(time, inits, paramsAll'); 
+   Ip = sim.variablevalues(:,1)'; 
+   Gp = sim.variablevalues(:,2)'; 
+   EGP = sim.variablevalues(:,5)';
+   S = sim.variablevalues(:,7)';
+   Ra = sim.variablevalues(:,11)';
+   U_id = sim.variablevalues(:,end)';
+   
+end
diff --git a/scripts/suport/IQMtools/IQMpro/tools/01-MEXmodels/CVODEMEX/lib/CVODEmex25.lib b/scripts/suport/IQMtools/IQMpro/tools/01-MEXmodels/CVODEMEX/lib/CVODEmex25.lib
index d0f44f34af0dfe7ca377c46f1bf1f8ab04af6927..d07a309494893288071813b86f6fe42674ae063a 100644
Binary files a/scripts/suport/IQMtools/IQMpro/tools/01-MEXmodels/CVODEMEX/lib/CVODEmex25.lib and b/scripts/suport/IQMtools/IQMpro/tools/01-MEXmodels/CVODEMEX/lib/CVODEmex25.lib differ