Merge pull request #10 from igem-waterloo/david/ane/0/ThermodynamicModel

David/ane/0/thermodynamic model

models/algae-n-enzymes/README.md

@@ -1 +1,30 @@
-# Algae 'n' Enzymes
+# Algae 'n' Enzymes
+
+## Summary
+For the Algae n' Enzymes project, we created a model to estimate the rate at which the C. vulgaris will be broken down in the cattle rumen, and thus estimated the amount of peir which will be released as a function of time and initial C. vulgaris fed to the cow. 
+A thermodynamics-based approach was used, based on [this paper](https://animres.edpsciences.org/articles/animres/abs/2006/05/z205012/z205012.html), assuming that the majority of the cell wall is composed of cellulose. Due to the limited and contradictory information regarding C. vulgaris cell wall concentrations in literature, a Monte Carlo algorithm was implemented for the glucose percentage of the cell wall. 
+
+### Data
+
+### Scripts
+
+MonteCarloMain.m: The main script containing the Monte Carlo model for glucose only.
+MultiMonteCarlo.m: A modification of the main script also allowing for the Monte Carlo simulation of time spent by feed in the rumen.
+odeModel.m: The ODE model from the paper.
+odeModelMulti.m: A modification of the ODE model to allow for varying time length inputs.
+
+### Assets
+The graphs resulting from a given set of runs if the Monte Carlo script were included. Note that these graphs need not necessarily coincide with future runs of the same script due to the random nature of the Monte Carlo simulation. 
+
+Average.jpg: 
+CellWallConcentration.jpg:
+DegradationOverTime.jpg:
+
+## Setup
+1. Make sure to install MATLAB and the ODE toolbox (which should be installed by default with your MATLAB installation) 
+2. Navigate to the models directory, by running 
+   `cd models\algae-n-enzymes\scripts`
+3. Ensure that the ODE model and the corresponding Monte Carlo script are located in the same directory on your local device 
+    (`MonteCarloMain.m` and `odeModel.m` for graphs or `MultiMonteCarlo.m` and `odeModelMulti.m` for numeric, time varying output)
+4. Run the model by running 
+   `MonteCarloMain.m` for the graphs or `MultiMonteCarlo.m` for the time-varying simulation.

models/algae-n-enzymes/scripts/MonteCarloMain.m

@@ -0,0 +1,111 @@
+clear all;
+clc;
+close all;
+
+numsimulations = 1000;
+resolution = 1000;
+AA = [zeros([numsimulations resolution])];
+BB = [zeros([numsimulations resolution])];
+T_sol = zeros(resolution);
+
+for i = 1:numsimulations
+
+% Define the number of random numbers
+n = 6;
+
+% Define the ranges for each number as [min, max]
+
+ranges = [
+    0.04, 0.54;  % Range for Rhamnose
+    0.004, 0.2;  % Range for Arabinose
+    0.01, 0.73; % Range for Glucose 
+    0, 0.26;  % Range for Galactose
+    0, 0.19; % Range for Xylose
+    0.025, 0.07   % Range for Mannose 
+];
+
+% Generate random numbers within the defined ranges
+randomNumbers = zeros(1, n);
+for k = 1:n
+    randomNumbers(k) = ranges(k, 1) + (ranges(k, 2) - ranges(k, 1)) * rand();
+end
+
+% Calculate the adjustment factor
+total = sum(randomNumbers);
+adjustmentFactor = 1 / total;
+
+% Adjust the numbers to sum to 1 while maintaining their proportion
+adjustedNumbers = randomNumbers * adjustmentFactor;
+
+RhaPer = adjustedNumbers(1);
+AraPer = adjustedNumbers(2);
+GluPer = adjustedNumbers(3);
+FucPer = adjustedNumbers(4);
+GalPer = adjustedNumbers(5);
+ManPer = adjustedNumbers(6);
+
+[t_sol,y_sol]=odeModel(GluPer, resolution);
+CW_sol = y_sol/y_sol(2);
+T_sol = t_sol;
+for q = 1:resolution
+AA(i,q)=CW_sol(q);
+BB(i,q)=y_sol(q);
+end
+end
+
+disp(AA);
+avg = mean(AA);
+
+
+figure(1)
+hold on
+title('Cell Wall Concentration')
+for l=1:size(BB,1)
+plot(T_sol, BB(l,:), 'r', 'DisplayName', 'CW_concentration')
+end
+hold off;
+
+figure(2)
+hold on
+title('Remaining Undegraded Cell Wall')
+for p=1:size(AA,1)
+plot(T_sol, AA(p,:), 'b')
+end
+ylabel('% Undegraded Cell Wall')
+xlabel('Time (min)')
+hold off;
+
+figure(3)
+title('Average Undegraded C. vulgaris Over Time')
+hold on;
+plot(T_sol, avg, 'g')
+ylabel('% Undegraded Cell Wall')
+xlabel('Time (min)')
+hold off;
+
+figure (4)
+hold on;
+title('Overlay')
+for p=1:size(AA,1)
+plot(T_sol, AA(p,:), 'b')
+end
+plot(T_sol, avg, 'g', 'LineWidth', 1.5)
+ylabel('% Undegraded Cell Wall')
+xlabel('Time (min)')
+hold off;
+
+figure (5)
+hold on;
+title('Degradation of C. vulgaris Over Time')
+for p=1:size(AA,1)
+plot(T_sol, 1-AA(p,:), 'b')
+end
+plot(T_sol, 1-avg, 'g', 'LineWidth', 1.5)
+ylabel('% Degraded C. vulgaris')
+xlabel('Time (min)')
+xlim([4 inf])
+h = zeros(2, 1);
+h(1) = plot(NaN,NaN,'b');
+h(2) = plot(NaN,NaN,'g');
+legend(h, 'Monte Carlo Runs', 'Average');
+hold off;

models/algae-n-enzymes/scripts/MultiMonteCarlo.m

@@ -0,0 +1,119 @@
+clear all;
+
+numsimulations = 50;
+resolution = 1000;
+AA = [zeros([numsimulations resolution])];
+BB = [zeros([numsimulations resolution])];
+T_sol = zeros([numsimulations resolution]);
+AV = [zeros([numsimulations 2])];
+
+for i = 1:numsimulations
+
+t_min = 60*24;
+t_max = 60*24*2;
+t_end = t_min+(t_max-t_min)*rand();
+%disp(t_end);
+% Define the number of random numbers
+n = 6;
+
+% Define the ranges for each number as [min, max]
+
+ranges = [
+    0.04, 0.54;  % Range for Rhamnose
+    0.004, 0.2;  % Range for Arabinose
+    0.01, 0.73; % Range for Glucose 
+    0, 0.26;  % Range for Galactose
+    0, 0.19; % Range for Xylose
+    0.025, 0.07   % Range for Mannose 
+];
+
+% Generate random numbers within the defined ranges
+randomNumbers = zeros(1, n);
+for k = 1:n
+    randomNumbers(k) = ranges(k, 1) + (ranges(k, 2) - ranges(k, 1)) * rand();
+end
+
+% Calculate the adjustment factor
+total = sum(randomNumbers);
+adjustmentFactor = 1 / total;
+
+% Adjust the numbers to sum to 1 while maintaining their proportion
+adjustedNumbers = randomNumbers * adjustmentFactor;
+
+RhaPer = adjustedNumbers(1);
+AraPer = adjustedNumbers(2);
+GluPer = adjustedNumbers(3);
+FucPer = adjustedNumbers(4);
+GalPer = adjustedNumbers(5);
+ManPer = adjustedNumbers(6);
+
+[t_sol,y_sol]=odeModelMulti(GluPer, resolution, t_end);
+%disp(size(t_sol));
+CW_sol = y_sol/y_sol(2);
+%disp(size(CW_sol'));
+%disp(size(AA))
+%plot(t_sol,CW_sol);
+for q = 1:resolution
+AA(i,q)=CW_sol(q);
+BB(i,q)=y_sol(q);
+T_sol(i,q) = t_sol(q);
+end
+AV(i,1) = t_sol(resolution);
+AV(i,2) = y_sol(resolution);
+end
+
+%disp(AA);
+avg = mean(AA);
+%disp(AV);
+disp("AVERAGE");
+disp(mean(AV));
+
+% figure(1)
+% hold on
+% title('Monte Carlo Concentration')
+% for l=1:size(BB,1)
+% plot(T_sol, BB(l,:), 'r', 'DisplayName', 'CW_concentration')
+% end
+% hold off;
+% 
+% figure(2)
+% hold on
+% title('Monte Carlo Percentage')
+% for p=1:size(AA,1)
+% plot(T_sol, AA(p,:), 'b')
+% end
+% ylabel('% Undegraded Cell Wall')
+% xlabel('Time (min)')
+% hold off;
+% 
+% figure(3)
+% title('Average')
+% hold on;
+% plot(T_sol, avg, 'g')
+% ylabel('% Undegraded Cell Wall')
+% xlabel('Time (min)')
+% hold off;
+% 
+% figure (4)
+% hold on;
+% title('Overlay')
+% for p=1:size(AA,1)
+% plot(T_sol(p,:), AA(p,:), 'b')
+% end
+% plot(T_sol, avg, 'g', 'LineWidth', 1.5)
+% ylabel('% Undegraded Cell Wall')
+% xlabel('Time (min)')
+% hold off;
+% 
+% figure (5)
+% hold on;
+% title('Degradation of C. vulgaris Over Time')
+% for p=1:size(AA,1)
+% plot(T_sol(p,:), 1-AA(p,:), 'b')
+% end
+% plot(T_sol, 1-avg, 'g', 'LineWidth', 1.5)
+% ylabel('% Degraded C. vulgaris')
+% xlabel('Time (min)')
+% xlim([5 1000])
+% % legend()
+% hold off;