diff --git a/AMF/+AMF/+utils/parfor_progress.m b/AMF/+AMF/+utils/parfor_progress.m new file mode 100644 index 0000000..07ec991 --- /dev/null +++ b/AMF/+AMF/+utils/parfor_progress.m @@ -0,0 +1,82 @@ +function percent = parfor_progress(N) +%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor. +% PARFOR_PROGRESS works by creating a file called parfor_progress.txt in +% your working directory, and then keeping track of the parfor loop's +% progress within that file. This workaround is necessary because parfor +% workers cannot communicate with one another so there is no simple way +% to know which iterations have finished and which haven't. +% +% PARFOR_PROGRESS(N) initializes the progress monitor for a set of N +% upcoming calculations. +% +% PARFOR_PROGRESS updates the progress inside your parfor loop and +% displays an updated progress bar. +% +% PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress +% bar. +% +% To suppress output from any of these functions, just ask for a return +% variable from the function calls, like PERCENT = PARFOR_PROGRESS which +% returns the percentage of completion. +% +% Example: +% +% N = 100; +% parfor_progress(N); +% parfor i=1:N +% pause(rand); % Replace with real code +% parfor_progress; +% end +% parfor_progress(0); +% +% See also PARFOR. + +% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/ + +% error(nargchk(0, 1, nargin, 'struct')); + +if nargin < 1 + N = -1; +end + +percent = 0; +w = 50; % Width of progress bar + +if N > 0 + f = fopen('parfor_progress.txt', 'w'); + if f<0 + error('Do you have write permissions for %s?', pwd); + end + fprintf(f, '%d\n', N); % Save N at the top of progress.txt + fclose(f); + + if nargout == 0 + disp([' 0%[>', repmat(' ', 1, w), ']']); + end +elseif N == 0 + delete('parfor_progress.txt'); + percent = 100; + + if nargout == 0 + disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']); + end +else + + if ~exist('parfor_progress.txt', 'file') + error('parfor_progress.txt not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize parfor_progress.txt.'); + end + + f = fopen('parfor_progress.txt', 'a'); + fprintf(f, '1\n'); + fclose(f); + + f = fopen('parfor_progress.txt', 'r'); + progress = fscanf(f, '%d'); + fclose(f); + percent = (length(progress)-1)/progress(1)*100; + + if nargout == 0 + perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage + disp([repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']); + end +end diff --git a/AMF/+AMF/+utils/parfor_progress2.m b/AMF/+AMF/+utils/parfor_progress2.m new file mode 100644 index 0000000..2335bb6 --- /dev/null +++ b/AMF/+AMF/+utils/parfor_progress2.m @@ -0,0 +1,82 @@ +function percent = parfor_progress2(N) +%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor. +% PARFOR_PROGRESS works by creating a file called parfor_progress.txt in +% your working directory, and then keeping track of the parfor loop's +% progress within that file. This workaround is necessary because parfor +% workers cannot communicate with one another so there is no simple way +% to know which iterations have finished and which haven't. +% +% PARFOR_PROGRESS(N) initializes the progress monitor for a set of N +% upcoming calculations. +% +% PARFOR_PROGRESS updates the progress inside your parfor loop and +% displays an updated progress bar. +% +% PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress +% bar. +% +% To suppress output from any of these functions, just ask for a return +% variable from the function calls, like PERCENT = PARFOR_PROGRESS which +% returns the percentage of completion. +% +% Example: +% +% N = 100; +% parfor_progress(N); +% parfor i=1:N +% pause(rand); % Replace with real code +% parfor_progress; +% end +% parfor_progress(0); +% +% See also PARFOR. + +% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/ + +% error(nargchk(0, 1, nargin, 'struct')); + +if nargin < 1 + N = -1; +end + +percent = 0; +w = 50; % Width of progress bar + +if N > 0 + f = fopen('parfor_progress.txt', 'w'); + if f<0 + error('Do you have write permissions for %s?', pwd); + end + fprintf(f, '%d\n', N); % Save N at the top of progress.txt + fclose(f); + + if nargout == 0 + disp([' 0%[>', repmat(' ', 1, w), ']']); + end +elseif N == 0 + delete('parfor_progress.txt'); + percent = 100; + + if nargout == 0 + disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']); + end +else + + if ~exist('parfor_progress.txt', 'file') + error('parfor_progress.txt not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize parfor_progress.txt.'); + end + +% f = fopen('parfor_progress.txt', 'a'); +% fprintf(f, '1\n'); +% fclose(f); + + f = fopen('parfor_progress.txt', 'r'); + progress = fscanf(f, '%d'); + fclose(f); + percent = (length(progress)-1)/progress(1)*100; + + if nargout == 0 + perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage + disp([perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']); + end +end diff --git a/AMF/+AMF/@ModelResult/plotAll.m b/AMF/+AMF/@ModelResult/plotAll.m index c2792a1..cbd97f3 100644 --- a/AMF/+AMF/@ModelResult/plotAll.m +++ b/AMF/+AMF/@ModelResult/plotAll.m @@ -1,6 +1,14 @@ -function this = plotAll(this, type, mode) +function this = plotAll(this, type, varargin) -if nargin < 3, mode = 'TRAJ'; end +if nargin < 3, mode = 'TRAJ'; +else mode = varargin{1}; +end + +if length(varargin)>1 + split = varargin{2}; +else + split = 9; +end import AMF.utils.defineCustomColormap @@ -9,15 +17,20 @@ comps = comps(logical([this.parameters.fit])); end + n = length(comps); -ns = sqrt(n); +ns = sqrt(split); +% ns = sqrt(n); numIter = this.options.numIter; figure('Name', upper(type)); for i = 1:n - subplot(ceil(ns),ceil(ns),i); hold on; + figure(ceil(i/split)); + counter = ceil(i/split); + subplot(ceil(ns),ceil(ns),i-(counter-1)*split); hold on; +% subplot(ceil(ns),ceil(ns),i); hold on; comp = comps(i); @@ -31,9 +44,10 @@ plotHist(this, comp, colorMap); case 'HIST_LOG' - - % TODO: plot logarithmic histograms - + % TODO: plot logarithmic histograms + case 'MAD' % median absolute deviation + plotMad(this, comp, 'g'); + otherwise error('Unknown plot mode %s', mode); end diff --git a/AMF/+AMF/@ModelResult/plotMad.m b/AMF/+AMF/@ModelResult/plotMad.m new file mode 100644 index 0000000..0bc42ce --- /dev/null +++ b/AMF/+AMF/@ModelResult/plotMad.m @@ -0,0 +1,22 @@ +function [hp1, hp2] = plotMad(this, comp, color) +% plot median absolute deviation +time_steps = this.options.numTimeSteps; +meta = zeros(time_steps,4); +for dt = 1:time_steps + y = quantile(comp.val(dt,:),[0.05 0.16 0.25 0.5 0.75 0.84 0.95], 2); + meta(dt,1) = y(4); + meta(dt,2) = median(abs(comp.val(dt,:)- meta(dt,1))); + meta(dt,3) = meta(dt,1)+meta(dt,2); + meta(dt,4) = meta(dt,1)-meta(dt,2); +end +time = 1:time_steps; +time = time*this.time(end)/time_steps; +X = [time fliplr(time)]; +Yupp = meta(:,3)'; +Ylwr = meta(:,4)'; +Y = [Ylwr fliplr(Yupp)]; +hp1 = fill(X, Y, color); +hp2 = plot(time, meta(:,1), '-k', 'LineWidth', 2); + +end + diff --git a/AMF/+AMF/@ModelResult/plot.m b/AMF/+AMF/@ModelResult/plots.m similarity index 56% rename from AMF/+AMF/@ModelResult/plot.m rename to AMF/+AMF/@ModelResult/plots.m index 49d5db2..8ded1a4 100644 --- a/AMF/+AMF/@ModelResult/plot.m +++ b/AMF/+AMF/@ModelResult/plots.m @@ -1,4 +1,4 @@ -function this = plot(this, names, mode) +function this = plots(this, names, mode) if nargin < 3, mode = 'TRAJ'; end @@ -31,6 +31,24 @@ case 'HIST_LOG' % TODO: plot logarithmic histograms + case 'MAD' % median absolute deviation + time_steps = this.options.numTimeSteps; + meta = zeros(time_steps,4); + for dt = 1:time_steps + y = quantile(comp.val(dt,:),[0.05 0.16 0.25 0.5 0.75 0.84 0.95], 2); + meta(dt,1) = y(4); + meta(dt,2) = median(abs(comp.val(dt,:)- meta(dt,1))); + meta(dt,3) = meta(dt,1)+meta(dt,2); + meta(dt,4) = meta(dt,1)-meta(dt,2); + end + time = 1:time_steps; + time = time*this.time(end)/time_steps; + X = [time fliplr(time)]; + Yupp = meta(:,3)'; + Ylwr = meta(:,4)'; + Y = [Ylwr fliplr(Yupp)]; + hp = fill(X, Y, 'g'); + plot(time, meta(:,1), '-k', 'LineWidth', 2); otherwise error('Unknown plot mode %s', mode); diff --git a/AMF/+AMF/runADAPT.m b/AMF/+AMF/runADAPT.m index 42224c8..7ca3d88 100644 --- a/AMF/+AMF/runADAPT.m +++ b/AMF/+AMF/runADAPT.m @@ -27,6 +27,7 @@ result = model.result; tic +utils.parfor_progress(numIter); parfor it = 1:numIter model2 = model; % parfor fix elt = 0; @@ -47,14 +48,17 @@ success = 1; catch err - % + disp(err.message); + utils.parfor_progress2; end end + utils.parfor_progress; - fprintf('Computed trajectory %d [%d] - %.2fs\n', it, max(model2.result.sse), elt); +% fprintf('Computed trajectory %d [%d] - %.2fs\n', it, max(model2.result.sse), elt); result(it) = model2.result; end +parfor_progress(0); toc result = AMF.ModelResult(model, result); diff --git a/AMF/models/dietModel/Hall/BW.emf b/AMF/models/dietModel/Hall/BW.emf new file mode 100644 index 0000000..5e27a8c Binary files /dev/null and b/AMF/models/dietModel/Hall/BW.emf differ diff --git a/AMF/models/dietModel/Hall/Hall.emf b/AMF/models/dietModel/Hall/Hall.emf new file mode 100644 index 0000000..88d01e0 Binary files /dev/null and b/AMF/models/dietModel/Hall/Hall.emf differ diff --git a/AMF/models/dietModel/Hall/Hall.m b/AMF/models/dietModel/Hall/Hall.m new file mode 100644 index 0000000..5d5f77d --- /dev/null +++ b/AMF/models/dietModel/Hall/Hall.m @@ -0,0 +1,9 @@ +function EE = Hall(t, dEI ) +%Energy Expenditure in mice according to Hall +LM = 25; +FM = -0.0649.*t.^2 + 2.0961.*t + 2.5263; + +EE = 2.1 + 0.331.*LM + 0.202.*FM + 0.4.*dEI; + +end + diff --git a/AMF/models/dietModel/Hall/Hall1.emf b/AMF/models/dietModel/Hall/Hall1.emf new file mode 100644 index 0000000..00418df Binary files /dev/null and b/AMF/models/dietModel/Hall/Hall1.emf differ diff --git a/AMF/models/dietModel/Hall/bw_fit.m b/AMF/models/dietModel/Hall/bw_fit.m new file mode 100644 index 0000000..7a0d333 --- /dev/null +++ b/AMF/models/dietModel/Hall/bw_fit.m @@ -0,0 +1,21 @@ +% Fit body weight data + +time = diets.t3/(24*7); % weeks HFD +TG = diets.perTG_m*853/1e6; % g TG + +% fit to polynomial +p = polyfit(time, TG, 2); % -0.0649*x^2 + 2.0961*x + 2.5263 + +% plot polynomial and data +x1 = linspace(0, 12, 100); +y1 = polyval(p,x1); + +figure(1);hold on; +plot(time, TG, 'ok', x1, y1, '-k', 'LineWidth', 2); +xlabel('Time on HFD (weeks)'); +ylabel('Body TG (g)'); +saveas(1, 'BW.emf'); + +% slope at t= 3 = 1.7067 g TG/week +% slope at t=10 = 0.7981 g TG/week + diff --git a/AMF/models/dietModel/Hall/hall_plot.m b/AMF/models/dietModel/Hall/hall_plot.m new file mode 100644 index 0000000..e8502a1 --- /dev/null +++ b/AMF/models/dietModel/Hall/hall_plot.m @@ -0,0 +1,11 @@ +% Hall +t = linspace(0,16,100); +t1 = [3 10]; +x1 = [11.3+1.35 12.9+1.35]; +ee3 = Hall(t, 3); +ee2 = Hall(t, 2); +ee4 = Hall(t, 4); +ee909 = Hall(t, 0.909); +plot(t,ee2, '--k', t, ee3, '-k', t, ee4, '--k',t, ee909, '--b', t1, x1, 'or', 'LineWidth', 2); +xlabel('Time on HFD (weeks)') +ylabel('Energy Expenditure (kcal/day)') \ No newline at end of file diff --git a/AMF/models/dietModel/Model_Diets_YP.m b/AMF/models/dietModel/Model_Diets_YP.m new file mode 100644 index 0000000..8954a88 --- /dev/null +++ b/AMF/models/dietModel/Model_Diets_YP.m @@ -0,0 +1,108 @@ +function MODEL = Model_Diets_YP() + +MODEL.DESCRIPTION = 'Model to explore glucose lipid interactions in the diet selection experiment.'; + +MODEL.PREDICTOR = { + 't' [0 24*7*12] {'time' 'hours' 'Time'} +}; + +MODEL.CONSTANTS = { +}; + +MODEL.PARAMETERS = { +% name fit init_value bnd meta-data + 'k1' 1 1 [] {} + 'k2' 1 1 [] {} + 'k3' 1 1 [] {} + 'k4' 1 1 [] {} + 'k5' 1 1 [] {} + 'k6' 1 1 [] {} + 'k7' 1 1 [] {} + 'k8' 1 1 [] {} + 'k9' 1 1 [] {} + 'k10' 1 1 [] {} + 'k11' 1 1 [] {} + 'k12' 1 1 [] {} + 'k13' 1 1 [] {} + 'k14' 1 1 [] {} + 'k15' 1 1 [] {} + 'k16' 1 1 [] {} + 'k17' 1 1 [] {} + 'k18' 1 1 [] {} + 'k19' 1 1 [] {} + 'k20' 1 1 [] {} + 'k21' 1 1 [] {} + 'k22' 1 1 [] {} + 'k23' 1 1 [] {} + 'k24' 1 1 [] {} + 'k25' 1 1 [] {} + 'k26' 1 1 [] {} + 'k27' 1 1 [] {} + 'k28' 1 1 [] {} + 'k29' 1 1 [] {} + 'k30' 1 1 [] {} + 'k31' 1 1 [] {} + 'k32' 1 1 [] {} +}; + +MODEL.STATES = { + 'hep_AcoA' 1 '3*v4 - v5 - v6' {} + 'hep_TG' 2 'v2 + v5/26 + v14/3 - v7 - v12 + v17' {} + 'hep_C' 1 'v6/13 + v15 + v16 - v8 - v11 - v13 + v32' {} + 'hep_BA' 0 'v8 - v9' {} + 'pl_FFA' 1 '3*v19 - v14' {} + 'pl_VLDL_TG' 1 'v12 - v20 - v31' {} + 'pl_VLDL_C' 1 'v13 - v15 - v21 - v22 + v18' {} + 'pl_HDL_C' 1 'v30 - v16 - v18' {} + 'per_AcoA' 1 '3*v26 - v27 - v28' {} + 'per_TG' 1 'v20 + v23 + v27/26 - v19 - v29' {} + 'per_C' 1 'v21 + v28/13 - v30' {} + 'pl_HDL_TG' 1 'v31 - v17' {} +}; + +MODEL.REACTIONS = { + % name reaction meta + 'v1' 'k1' {} + 'v2' 'k2' {} + 'v3' 'k3' {} + 'v4' 'k4' {} + 'v5' 'k5 * hep_AcoA' {} + 'v6' 'k6 * hep_AcoA' {} + 'v7' 'k7 * hep_TG' {} + 'v8' 'k8 * hep_C' {} + 'v9' 'k9 * hep_BA' {} + 'v10' 'k10 * hep_BA' {} + 'v11' 'k11 * hep_C' {} + 'v12' 'k12 * hep_TG' {} + 'v13' 'k13 * hep_C' {} + 'v14' 'k14 * pl_FFA' {} + 'v15' 'k15 * pl_VLDL_C' {} + 'v16' 'k16 * pl_HDL_C' {} + 'v17' 'k17 * pl_HDL_TG' {} + 'v18' 'k18 * (pl_HDL_C / (pl_HDL_TG + pl_HDL_C) - pl_VLDL_C / (pl_VLDL_C + pl_VLDL_TG))' {} % HDL_C -> VLDL_C + 'v19' 'k19 * per_TG' {} + 'v20' 'k20 * pl_VLDL_TG' {} + 'v21' 'k21 * pl_VLDL_C' {} + 'v22' 'k22 * pl_VLDL_C' {} + 'v23' 'k23' {} + 'v24' 'k24' {} + 'v25' 'k25' {} + 'v26' 'k26' {} + 'v27' 'k27 * per_AcoA' {} + 'v28' 'k28 * per_AcoA' {} + 'v29' 'k29 * per_TG' {} + 'v30' 'k30 * per_C' {} + 'v31' 'k31 * (pl_VLDL_TG / (pl_VLDL_TG + pl_VLDL_C) - pl_HDL_TG / (pl_HDL_TG + pl_HDL_C))' {} % VLDL_TG -> HDL_TG + 'v32' 'k32' {} % cholesterol from diet to liver + + % helper functions + 'pl_TC' 'pl_VLDL_C + pl_HDL_C' {} + 'fat_intk' 'v1 + v2 + v23 + v24' {} + 'gluc_intk' 'v3 + v4 + v25 + v26' {} + 'fat_ox' 'v1 + v24 + v7 + v29' {} + 'gluc_ox' 'v3 + v25' {} +% 'chol_exc' 'v42 - v41 + v11' {} + + +}; + diff --git a/AMF/models/dietModel/dietData.m b/AMF/models/dietModel/dietData.m new file mode 100644 index 0000000..96efbfa --- /dev/null +++ b/AMF/models/dietModel/dietData.m @@ -0,0 +1,30 @@ +function DATASET = dietData() + +DATASET.DESCRIPTION = 'Diet data.'; + +DATASET.FILE = 'dietData'; + +DATASET.GROUPS = { + 'diets' +% 'diets_R' +% 'diets_NAMR' +}; + +DATASET.FIELDS = { + %name %obs %t %mean % stand dev % unit conv + 'hep_TG' 1 't1' 'hep_TG_m' 'hep_TG_std' 1 [] + 'hep_C' 1 't1' 'hep_C_m' 'hep_C_std' 1 [] + 'pl_VLDL_TG' 1 't2' 'pl_VLDL_TG_m' 'pl_VLDL_TG_std' 1 [] + 'pl_HDL_C' 1 't2' 'pl_HDL_C_m' 'pl_HDL_C_std' 1 [] + 'pl_FFA' 1 't2' 'pl_FFA_m' 'pl_FFA_std' 1 [] + 'per_TG' 1 't3' 'perTG_m' 'perTG_std' 1 [] + 'pl_TC' 1 't2' 'pl_TC_m' 'pl_TC_std' 1 [] + 'fat_intk' 1 't2' 'FatInt_m' 'FatInt_std' 1 [] + 'gluc_intk' 1 't2' 'GlucInt_m' 'GlucInt_std' 1 [] + 'fat_ox' 1 't4' 'fat_ox_m' 'fat_ox_std' 1 [] + 'gluc_ox' 1 't4' 'gluc_ox_m' 'gluc_ox_std' 1 [] +% 'chol_intk' 0 't2' 'chol_intk_m' 'chol_intk_std' 1 [] +}; + +DATASET.FUNCTIONS = { +}; \ No newline at end of file diff --git a/AMF/models/dietModel/dietData.mat b/AMF/models/dietModel/dietData.mat new file mode 100644 index 0000000..2420587 Binary files /dev/null and b/AMF/models/dietModel/dietData.mat differ diff --git a/AMF/models/dietModel/loadData.m b/AMF/models/dietModel/loadData.m new file mode 100644 index 0000000..7d2dddb --- /dev/null +++ b/AMF/models/dietModel/loadData.m @@ -0,0 +1,153 @@ +data.diets.t1 = [0 12]*7*24; +data.diets.t2 = [0 4 8 12]*7*24; +data.diets.t3 = [0 1 3 5 7 9 12]*7*24; +data.diets.t4 = [3 10]*7*24; + + +% liver weights (g) +liver_weights = [ 2.85 + 3.40 + 3.39 + 1.88 + 1.74 + 1.43 + 3.02 + 3.96 + 4.28 + 3.54 ]; + +% nmol/mg tissue (mumol/g tissue) +data.diets.hep_TG_r = [ 109.3 + 124.6 + 115.9 + 62.5 + 44.0 + 47.9 + 106.0 + 111.9 + 112.1 + 83.8 ]; +data.diets.hep_TG_m = nanmean(data.diets.hep_TG_r .* liver_weights); +data.diets.hep_TG_std = nanstd(data.diets.hep_TG_r .* liver_weights); +% Wielinga et al. 9.1 (0.8) mumol TG/liver +data.diets.hep_TG_m = [9.1 data.diets.hep_TG_m]; % mumol / liver +data.diets.hep_TG_std = [0.8 data.diets.hep_TG_std]; % mumol / liver +% data.diets.hep_TG_m = [8.1 91.79]; + +% nmol/mg tissue (mumol/g tissue) +data.diets.hep_C_r = [ 16.6 + 12.3 + 14.7 + 8.8 + 9.5 + 10.2 + 13.5 + 15.7 + 13.3 + 13.2 ]; +data.diets.hep_C_m = nanmean(data.diets.hep_C_r .* liver_weights); +data.diets.hep_C_std = nanstd(data.diets.hep_C_r .* liver_weights); +% Wielinga et al. 6.8 (0.4) mumol C/liver +data.diets.hep_C_m = [6.8 data.diets.hep_C_m]; % mumol / liver +data.diets.hep_C_std = [0.4 data.diets.hep_C_std]; % mumol / liver +% data.diets.hep_C_m = [6.1 12.79]; +% data.diets.hep_C_std = [0.35 2.61]; + +% mM (mumol, assuming plasma volume = 1ml) +data.diets.pl_FFA_r = [ 1.87 0.71 0.53 1.40 + 1.27 0.87 0.55 1.51 + 1.52 0.65 0.61 1.64 + 0.99 0.66 0.65 1.04 + 0.77 0.61 0.46 1.19 + 1.60 0.61 0.38 1.07 + 1.34 0.52 0.80 1.13 + 1.41 0.92 0.59 1.17 + 1.51 0.68 0.46 1.53 + 1.43 0.61 0.56 1.48 ]; +data.diets.pl_FFA_m = [ 1.37 0.68 0.56 1.32]; % mumol / plasma +data.diets.pl_FFA_std = [0.31 0.12 0.12 0.22]; % mumol / plasma + +% mM +data.diets.pl_VLDL_TG_r = [ 2.27 1.58 1.45 2.00 + 2.7 1.25 1.47 2.39 + 2.78 1.08 1.31 2.42 + 2.07 1.30 1.10 0.77 + 1.88 1.06 1.17 0.79 + 5.72 1.64 1.13 0.82 + 4.23 1.57 1.34 2.16 + 1.89 1.20 2.30 7.97 + 2.57 1.44 2.48 7.50 + 1.89 1.39 2.02 5.57] ; + +data.diets.pl_VLDL_TG_m = [2.80 1.35 1.58 3.24]; % mumol / plasma +data.diets.pl_VLDL_TG_std = [1.25 0.21 0.50 2.75]; % mumol / plasma + +% mM +data.diets.pl_TC_r = [ 3.17 5.01 4.94 5.79 + 3.22 4.39 4.92 7.07 + 3.67 4.32 4.94 6.86 + 2.45 4.22 4.26 5.25 + 2.20 3.92 3.96 5.13 + 5.25 4.95 3.37 4.52 + 4.56 4.74 4.64 7.09 + 2.82 5.31 6.21 15.45 + 3.07 5.31 6.06 14.58 + 2.55 4.89 5.44 15.22] ; +data.diets.pl_TC_m = [ 3.30 4.70 4.87 8.69]; % mumol / plasma +data.diets.pl_TC_std = [0.96 0.47 0.88 4.50]; % mumol / plasma + +% mM +data.diets.pl_HDL_C_r = [ 1.69 3.03 3.21 2.64 + 1.29 3.03 3.44 2.47 + 1.06 3.35 3.25 2.51 + 0.61 3.09 3.64 3.56 + 0.82 2.69 2.83 2.73 + 1.79 1.82 NaN 2.86 + 1.10 1.91 3.91 2.97 + 1.31 2.56 2.63 1.41 + 1.15 3.03 2.86 1.81 + 1.20 NaN 3.49 1.80]; +data.diets.pl_HDL_C_m = [1.20 2.72 3.25 2.48]; % mumol / plasma +data.diets.pl_HDL_C_std = [0.36 0.54 0.42 0.64]; % mumol / plasma + + +% TG = (fat mass (g)) / 853 (g/mol) = TG (mol) ; + +data.diets.fat_mass = [ 2.4 7.6 11.3 14.4 14.4 18.1 19.0 + 3.6 8.1 10.2 13.8 15.7 18.4 18.6 + 2.3 6.3 10.1 13.7 15.8 19.4 20.8 + 0.8 2.9 5.0 8.4 10.0 13.4 13.8 + 0.9 2.8 3.9 7.2 8.4 10.7 11.7 + 1.1 1.8 1.6 3.4 4.5 6.7 9.6 + 1.3 3.2 4.7 9.3 11.9 15.3 19.9 + 3.2 8.2 11.8 14.7 16.7 21.1 23.3 + 2.1 6.9 10.8 14.7 17.5 21.3 22.4 + 1.7 6.2 11.6 15.4 18.3 22.4 23.4]; + +data.diets.perTG_m = nanmean(data.diets.fat_mass*(1e6/853)); +data.diets.perTG_std = nanstd(data.diets.fat_mass*(1e6/853)); + +% energy intake in TG equivalents +% TG = 9.2 kcal/g +% fat intake +data.diets.FatInt_m = [3.1 3.1 3.1 3.1]*0.60*(5.2/9.2)*1*(1e6/853)/24; +data.diets.FatInt_std = [0.3 0.3 0.3 0.3]*0.60*(5.2/9.2)*1*(1e6/853)/24; + +% gluc = 4.2 kcal/g --> gluc / TG = 0.46; 4/9 = 0.44; +% glucose intake +data.diets.GlucInt_m = [3.1 3.1 3.1 3.1]*0.20*(5.2/4.2)*1*(1e6/180)/24; +data.diets.GlucInt_std = [0.3 0.3 0.3 0.3]*0.20*(5.2/4.2)*1*(1e6/180)/24; + +% chol intake %% is there a TG equivalent of cholesterol? +data.diets.chol_intk_m = [3.1 3.1 3.1 3.1]*0.01; +data.diets.chol_intk_std = [0.3 0.3 0.3 0.3]*0.01; + +% fat oxidation (kcal/h) +data.diets.fat_ox_m = [6.98 8.79]*(1e6/853)/9.2/24; +data.diets.fat_ox_std = [0.81 0.78]*(1e6/853)/9.2/24; + +% glucose oxidation (kcal/h) +data.diets.gluc_ox_m = [1.32 1.10]*(1e6/180)/4.2/24; +data.diets.gluc_ox_std = [0.19 0.09]*(1e6/180)/4.2/24; + +save('dietData.mat', '-struct', 'data'); \ No newline at end of file diff --git a/AMF/models/dietModel/results/_Model_Diets_YP_diets__10_50.mat b/AMF/models/dietModel/results/_Model_Diets_YP_diets__10_50.mat new file mode 100644 index 0000000..6503f86 Binary files /dev/null and b/AMF/models/dietModel/results/_Model_Diets_YP_diets__10_50.mat differ diff --git a/AMF/models/dietModel/runDiet.m b/AMF/models/dietModel/runDiet.m new file mode 100644 index 0000000..50245b0 --- /dev/null +++ b/AMF/models/dietModel/runDiet.m @@ -0,0 +1,37 @@ +%% initialize + +import AMF.* + +model = Model('Model_Diets_YP'); +data = DataSet('dietData'); + +loadGroup(data, 'diets'); +initiateExperiment(model, data); + +%% config + +model.options.useMex = 1; +model.options.savePrefix = ''; +model.options.odeTol = [1e-12 1e-12 100]; +model.options.numIter = 10; +model.options.numTimeSteps = 50; +model.options.parScale = [2 -2]; +model.options.seed = 1; +model.options.SSTime = 1000; +model.options.lab1 = .1; + +parseAll(model); +compileAll(model); + +%% run + +result = runADAPT(model); + + +%% plot + +% plotAll(result, 'parameters', 'traj'); +plotAll(result, 'states', 'mad'); +plotAll(result, 'reactions', 'mad'); + +% random change \ No newline at end of file diff --git a/AMF/models/dietModel/temp/C_Model_Diets_YP_ODE.m b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODE.m new file mode 100644 index 0000000..1cabf2f --- /dev/null +++ b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODE.m @@ -0,0 +1,102 @@ +function dxdt = C_Model_Diets_YP_ODE(t,x,p,u,m) + + + +hep_AcoA = x(m.s.hep_AcoA); +hep_TG = x(m.s.hep_TG); +hep_C = x(m.s.hep_C); +hep_BA = x(m.s.hep_BA); +pl_FFA = x(m.s.pl_FFA); +pl_VLDL_TG = x(m.s.pl_VLDL_TG); +pl_VLDL_C = x(m.s.pl_VLDL_C); +pl_HDL_C = x(m.s.pl_HDL_C); +per_AcoA = x(m.s.per_AcoA); +per_TG = x(m.s.per_TG); +per_C = x(m.s.per_C); +pl_HDL_TG = x(m.s.pl_HDL_TG); + +k1 = p(m.p.k1); +k2 = p(m.p.k2); +k3 = p(m.p.k3); +k4 = p(m.p.k4); +k5 = p(m.p.k5); +k6 = p(m.p.k6); +k7 = p(m.p.k7); +k8 = p(m.p.k8); +k9 = p(m.p.k9); +k10 = p(m.p.k10); +k11 = p(m.p.k11); +k12 = p(m.p.k12); +k13 = p(m.p.k13); +k14 = p(m.p.k14); +k15 = p(m.p.k15); +k16 = p(m.p.k16); +k17 = p(m.p.k17); +k18 = p(m.p.k18); +k19 = p(m.p.k19); +k20 = p(m.p.k20); +k21 = p(m.p.k21); +k22 = p(m.p.k22); +k23 = p(m.p.k23); +k24 = p(m.p.k24); +k25 = p(m.p.k25); +k26 = p(m.p.k26); +k27 = p(m.p.k27); +k28 = p(m.p.k28); +k29 = p(m.p.k29); +k30 = p(m.p.k30); +k31 = p(m.p.k31); +k32 = p(m.p.k32); + +v1 = k1; +v2 = k2; +v3 = k3; +v4 = k4; +v5 = k5 * hep_AcoA; +v6 = k6 * hep_AcoA; +v7 = k7 * hep_TG; +v8 = k8 * hep_C; +v9 = k9 * hep_BA; +v10 = k10 * hep_BA; +v11 = k11 * hep_C; +v12 = k12 * hep_TG; +v13 = k13 * hep_C; +v14 = k14 * pl_FFA; +v15 = k15 * pl_VLDL_C; +v16 = k16 * pl_HDL_C; +v17 = k17 * pl_HDL_TG; +v18 = k18 * (pl_HDL_C / (pl_HDL_TG + pl_HDL_C) - pl_VLDL_C / (pl_VLDL_C + pl_VLDL_TG)); +v19 = k19 * per_TG; +v20 = k20 * pl_VLDL_TG; +v21 = k21 * pl_VLDL_C; +v22 = k22 * pl_VLDL_C; +v23 = k23; +v24 = k24; +v25 = k25; +v26 = k26; +v27 = k27 * per_AcoA; +v28 = k28 * per_AcoA; +v29 = k29 * per_TG; +v30 = k30 * per_C; +v31 = k31 * (pl_VLDL_TG / (pl_VLDL_TG + pl_VLDL_C) - pl_HDL_TG / (pl_HDL_TG + pl_HDL_C)); +v32 = k32; +pl_TC = pl_VLDL_C + pl_HDL_C; +fat_intk = v1 + v2 + v23 + v24; +gluc_intk = v3 + v4 + v25 + v26; +fat_ox = v1 + v24 + v7 + v29; +gluc_ox = v3 + v25; + +dxdt(1) = 3*v4 - v5 - v6; +dxdt(2) = v2 + v5/26 + v14/3 - v7 - v12 + v17; +dxdt(3) = v6/13 + v15 + v16 - v8 - v11 - v13 + v32; +dxdt(4) = v8 - v9; +dxdt(5) = 3*v19 - v14; +dxdt(6) = v12 - v20 - v31; +dxdt(7) = v13 - v15 - v21 - v22 + v18; +dxdt(8) = v30 - v16 - v18; +dxdt(9) = 3*v26 - v27 - v28; +dxdt(10) = v20 + v23 + v27/26 - v19 - v29; +dxdt(11) = v21 + v28/13 - v30; +dxdt(12) = v31 - v17; + +dxdt = dxdt(:); \ No newline at end of file diff --git a/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEC.mexw64 b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEC.mexw64 new file mode 100644 index 0000000..a613c8c Binary files /dev/null and b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEC.mexw64 differ diff --git a/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEMEX.m b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEMEX.m new file mode 100644 index 0000000..31b93e2 --- /dev/null +++ b/AMF/models/dietModel/temp/C_Model_Diets_YP_ODEMEX.m @@ -0,0 +1,102 @@ +function dxdt = C_Model_Diets_YP_ODEMEX(t,x,p,u,m) + + + +hep_AcoA = x(m.s.hep_AcoA); +hep_TG = x(m.s.hep_TG); +hep_C = x(m.s.hep_C); +hep_BA = x(m.s.hep_BA); +pl_FFA = x(m.s.pl_FFA); +pl_VLDL_TG = x(m.s.pl_VLDL_TG); +pl_VLDL_C = x(m.s.pl_VLDL_C); +pl_HDL_C = x(m.s.pl_HDL_C); +per_AcoA = x(m.s.per_AcoA); +per_TG = x(m.s.per_TG); +per_C = x(m.s.per_C); +pl_HDL_TG = x(m.s.pl_HDL_TG); + +k1 = p(m.p.k1); +k2 = p(m.p.k2); +k3 = p(m.p.k3); +k4 = p(m.p.k4); +k5 = p(m.p.k5); +k6 = p(m.p.k6); +k7 = p(m.p.k7); +k8 = p(m.p.k8); +k9 = p(m.p.k9); +k10 = p(m.p.k10); +k11 = p(m.p.k11); +k12 = p(m.p.k12); +k13 = p(m.p.k13); +k14 = p(m.p.k14); +k15 = p(m.p.k15); +k16 = p(m.p.k16); +k17 = p(m.p.k17); +k18 = p(m.p.k18); +k19 = p(m.p.k19); +k20 = p(m.p.k20); +k21 = p(m.p.k21); +k22 = p(m.p.k22); +k23 = p(m.p.k23); +k24 = p(m.p.k24); +k25 = p(m.p.k25); +k26 = p(m.p.k26); +k27 = p(m.p.k27); +k28 = p(m.p.k28); +k29 = p(m.p.k29); +k30 = p(m.p.k30); +k31 = p(m.p.k31); +k32 = p(m.p.k32); + +v1 = k1; +v2 = k2; +v3 = k3; +v4 = k4; +v5 = k5 * hep_AcoA; +v6 = k6 * hep_AcoA; +v7 = k7 * hep_TG; +v8 = k8 * hep_C; +v9 = k9 * hep_BA; +v10 = k10 * hep_BA; +v11 = k11 * hep_C; +v12 = k12 * hep_TG; +v13 = k13 * hep_C; +v14 = k14 * pl_FFA; +v15 = k15 * pl_VLDL_C; +v16 = k16 * pl_HDL_C; +v17 = k17 * pl_HDL_TG; +v18 = k18 * (pl_HDL_C / (pl_HDL_TG + pl_HDL_C) - pl_VLDL_C / (pl_VLDL_C + pl_VLDL_TG)); +v19 = k19 * per_TG; +v20 = k20 * pl_VLDL_TG; +v21 = k21 * pl_VLDL_C; +v22 = k22 * pl_VLDL_C; +v23 = k23; +v24 = k24; +v25 = k25; +v26 = k26; +v27 = k27 * per_AcoA; +v28 = k28 * per_AcoA; +v29 = k29 * per_TG; +v30 = k30 * per_C; +v31 = k31 * (pl_VLDL_TG / (pl_VLDL_TG + pl_VLDL_C) - pl_HDL_TG / (pl_HDL_TG + pl_HDL_C)); +v32 = k32; +pl_TC = pl_VLDL_C + pl_HDL_C; +fat_intk = v1 + v2 + v23 + v24; +gluc_intk = v3 + v4 + v25 + v26; +fat_ox = v1 + v24 + v7 + v29; +gluc_ox = v3 + v25; + +dxdt(1) = 3*v4 - v5 - v6; +dxdt(2) = v2 + v5/26 + v14/3 - v7 - v12 + v17; +dxdt(3) = v6/13 + v15 + v16 - v8 - v11 - v13 + v32; +dxdt(4) = v8 - v9; +dxdt(5) = 3*v19 - v14; +dxdt(6) = v12 - v20 - v31; +dxdt(7) = v13 - v15 - v21 - v22 + v18; +dxdt(8) = v30 - v16 - v18; +dxdt(9) = 3*v26 - v27 - v28; +dxdt(10) = v20 + v23 + v27/26 - v19 - v29; +dxdt(11) = v21 + v28/13 - v30; +dxdt(12) = v31 - v17; + +dxdt = dxdt(:); \ No newline at end of file diff --git a/AMF/models/dietModel/temp/C_Model_Diets_YP_REACTIONS.m b/AMF/models/dietModel/temp/C_Model_Diets_YP_REACTIONS.m new file mode 100644 index 0000000..37bdc9b --- /dev/null +++ b/AMF/models/dietModel/temp/C_Model_Diets_YP_REACTIONS.m @@ -0,0 +1,125 @@ +function v = C_Model_Diets_YP_REACTIONS(t,x,p,u,m) + + + +hep_AcoA = x(m.s.hep_AcoA); +hep_TG = x(m.s.hep_TG); +hep_C = x(m.s.hep_C); +hep_BA = x(m.s.hep_BA); +pl_FFA = x(m.s.pl_FFA); +pl_VLDL_TG = x(m.s.pl_VLDL_TG); +pl_VLDL_C = x(m.s.pl_VLDL_C); +pl_HDL_C = x(m.s.pl_HDL_C); +per_AcoA = x(m.s.per_AcoA); +per_TG = x(m.s.per_TG); +per_C = x(m.s.per_C); +pl_HDL_TG = x(m.s.pl_HDL_TG); + +k1 = p(m.p.k1); +k2 = p(m.p.k2); +k3 = p(m.p.k3); +k4 = p(m.p.k4); +k5 = p(m.p.k5); +k6 = p(m.p.k6); +k7 = p(m.p.k7); +k8 = p(m.p.k8); +k9 = p(m.p.k9); +k10 = p(m.p.k10); +k11 = p(m.p.k11); +k12 = p(m.p.k12); +k13 = p(m.p.k13); +k14 = p(m.p.k14); +k15 = p(m.p.k15); +k16 = p(m.p.k16); +k17 = p(m.p.k17); +k18 = p(m.p.k18); +k19 = p(m.p.k19); +k20 = p(m.p.k20); +k21 = p(m.p.k21); +k22 = p(m.p.k22); +k23 = p(m.p.k23); +k24 = p(m.p.k24); +k25 = p(m.p.k25); +k26 = p(m.p.k26); +k27 = p(m.p.k27); +k28 = p(m.p.k28); +k29 = p(m.p.k29); +k30 = p(m.p.k30); +k31 = p(m.p.k31); +k32 = p(m.p.k32); + +v1 = k1; +v2 = k2; +v3 = k3; +v4 = k4; +v5 = k5 * hep_AcoA; +v6 = k6 * hep_AcoA; +v7 = k7 * hep_TG; +v8 = k8 * hep_C; +v9 = k9 * hep_BA; +v10 = k10 * hep_BA; +v11 = k11 * hep_C; +v12 = k12 * hep_TG; +v13 = k13 * hep_C; +v14 = k14 * pl_FFA; +v15 = k15 * pl_VLDL_C; +v16 = k16 * pl_HDL_C; +v17 = k17 * pl_HDL_TG; +v18 = k18 * (pl_HDL_C / (pl_HDL_TG + pl_HDL_C) - pl_VLDL_C / (pl_VLDL_C + pl_VLDL_TG)); +v19 = k19 * per_TG; +v20 = k20 * pl_VLDL_TG; +v21 = k21 * pl_VLDL_C; +v22 = k22 * pl_VLDL_C; +v23 = k23; +v24 = k24; +v25 = k25; +v26 = k26; +v27 = k27 * per_AcoA; +v28 = k28 * per_AcoA; +v29 = k29 * per_TG; +v30 = k30 * per_C; +v31 = k31 * (pl_VLDL_TG / (pl_VLDL_TG + pl_VLDL_C) - pl_HDL_TG / (pl_HDL_TG + pl_HDL_C)); +v32 = k32; +pl_TC = pl_VLDL_C + pl_HDL_C; +fat_intk = v1 + v2 + v23 + v24; +gluc_intk = v3 + v4 + v25 + v26; +fat_ox = v1 + v24 + v7 + v29; +gluc_ox = v3 + v25; + +v(1) = v1; +v(2) = v2; +v(3) = v3; +v(4) = v4; +v(5) = v5; +v(6) = v6; +v(7) = v7; +v(8) = v8; +v(9) = v9; +v(10) = v10; +v(11) = v11; +v(12) = v12; +v(13) = v13; +v(14) = v14; +v(15) = v15; +v(16) = v16; +v(17) = v17; +v(18) = v18; +v(19) = v19; +v(20) = v20; +v(21) = v21; +v(22) = v22; +v(23) = v23; +v(24) = v24; +v(25) = v25; +v(26) = v26; +v(27) = v27; +v(28) = v28; +v(29) = v29; +v(30) = v30; +v(31) = v31; +v(32) = v32; +v(33) = pl_TC; +v(34) = fat_intk; +v(35) = gluc_intk; +v(36) = fat_ox; +v(37) = gluc_ox; diff --git a/AMF/models/minGlucModel2/runMinGluc2.m b/AMF/models/minGlucModel2/runMinGluc2.m index 999e3ca..9c9898f 100644 --- a/AMF/models/minGlucModel2/runMinGluc2.m +++ b/AMF/models/minGlucModel2/runMinGluc2.m @@ -11,10 +11,11 @@ %% config -model.options.optimset.Display = 'off'; -model.options.useMex = 1; model.options.numIter = 500; model.options.numTimeSteps = 100; + +model.options.optimset.Display = 'off'; +model.options.useMex = 1; model.options.SSTime = 30; model.options.lab1 = .1; model.options.randPars = 0; @@ -25,10 +26,11 @@ %% run -% model.functions.reg = @minGlucReg; - result = runADAPT(model); + +%check whether AUC satisfies carbohydrate dose + %% plot close all