Difference between revisions of "Team:Dundee/Modeling/Appendix1"
Fmacfarlane (Talk | contribs) |
Fmacfarlane (Talk | contribs) |
||
| Line 401: | Line 401: | ||
</section> | </section> | ||
| − | <div class="content"><pre class="codeinput"><span class="comment">%Function to define the non-dimensionalised system of ODEs describing the</span> | + | <div class="content"><pre class="codeinput"><span class="comment"><h2>MATLAB code from obp.m file</h2>%Function to define the non-dimensionalised system of ODEs describing the</span> |
<span class="comment">%folding of the Oderant Binding Protein (OBP). u is a 3 dimensional vector</span> | <span class="comment">%folding of the Oderant Binding Protein (OBP). u is a 3 dimensional vector</span> | ||
<span class="comment">%where;</span> | <span class="comment">%where;</span> | ||
Revision as of 09:47, 10 August 2015
Blood: Haptoglobin and Haemoglobin Binding MATLAB Code
Before using MATLAB the original system of equations was non-dimensionalised to the system:
$$ \begin{eqnarray*} \frac{dHp}{d\tau}&=&[Hp\cdot \alpha_{H}]-\lambda Hp \alpha_{H},\\ \frac{d\alpha_{H}}{d\tau}&=&[Hp\cdot \alpha_{H}]-\lambda Hp \alpha_{H},\\ \frac{d[Hp \cdot \alpha_{H}]}{d\tau}&=&\lambda Hp \alpha_{H}-[Hp \cdot \alpha_{H}]-\gamma [Hp \cdot \alpha_{H}],\\ \frac{d[Hp \cdot \alpha_{H} \cdot \beta_{H}]}{d\tau}&=&\gamma [Hp \cdot \alpha_{H}]. \end{eqnarray*} $$The initial conditions were also non-dimensionalised to become:
$$ \begin{eqnarray*} Hp(0)&=&4.17,\\ \alpha_{H}(0)&=&1,\\ [Hp \cdot \alpha_{H}](0)&=&0,\\ [Hp \cdot \alpha_{H} \cdot \beta_{H}](0)&=&0. \end{eqnarray*} $$To perform senstivity analysis that is described in the Blood section of the BioSpray models, MATLAB was used. Two files were written one to set the function, haptosen.m, and one to solve the function and plot the results in Figure 1, run_haptosen.m. Both files are shown below, where the green is comments to aid in understanding of the scripts.
MATLAB code from haptosen.m file
%Function to define the non-dimensionalised system of ODEs describing the %binding between Haptoglobin and Haemoglobin. u is a 4 dimensional vector %where; % u(1)=Hp (Haptoglobin concentration) % u(2)=alphaHemo (Free Hemoglobin) % u(3)=alphaHemo/Hapto complex (Haptoglobin + alphaHemo complex) % u(4)=Hemo/Hapto complex (Full hemo/hapto complex). % t is the time that the simulation is run over. % lambda represents the parameter in the system determined by binding rates nd initial concentration of haemoglobin. % gamma represents the parameter in the system determined by binding rates. function f = haptosen(t,u,lambda,gamma); % define the system of ODEs by setting the vector f = left hand side of the system. f = [u(3)-lambda.*u(1)*u(2); u(3)-lambda.*u(1)*u(2); -u(3)+lambda.*u(1)*u(2)-gamma.*u(3); gamma.*u(3)]; % The function is then ended. end % This function can then be called in the run_haptosen.m file by using % @haptosen command.MATLAB code from run_haptosen.m file
File to perform sensitivity analysis for Haptoglobin and Haemoglobin % binding model. % a is the number of values chosen for each of the ranges of parameters. % Therefore 100 different values were chosen. a=100; %Lambda1 is the range of values for the parameter Lambda, chosen with the %expected value of (100/317)*0.3878494524 as the mean value. lambda1=linspace(0,((100/317)*0.3878494524)*2,a); %Gamma1 is the range of values for the parameter Gamma, chosen with the %expected value of (83/317) as the mean value. gamma1=linspace(0,(83/317)*2,a); %T is the final time of the simulation. T=30; % Store sets an empty matrix for the final concentrations of the complex % with the varied parameter values. Store = zeros(a,a); %figure(1) calls up a new figure. figure(1) % The for loop solves the ode system defined in haptosen function using ode23 for the range of values for both % parameters. for i=1:a for j=1:a [t,u]=ode23(@haptosen,[0 T],[4.17 1 0 0],[],lambda1(a+1-i),gamma1(j)); Store(i,j) = u(end,4); % Store(i,j) = (u(3,4)-u(2,4))/(t(3)-t(2)); end end %imagesc plots the surf plot of the parameters and the complex %concentration. imagesc(Store) % xlabel and ylabel add labels to the plot. xlabel('Increasing \lambda','FontSize',15,'FontWeight','bold'); ylabel('Increasing \gamma','FontSize',15,'FontWeight','bold'); %set removes the x and y axis set(gca,'YTick',[],'XTick',[]); % The following lines add a colour bar with defined labal and text. c=colorbar('Ticks',[0,0.99],'TickLabels',{'None','High'},'FontSize',15,'FontWeight','bold'); c.Label.String = 'Complex Concentration'; %the following lines add a text arrow to allow for annotation of the graph with defined colours positions and width. ta1 = annotation('textarrow', [0.13 0.79], [0.13 0.92]); h = text(0.5,0.5,'Complex Formation'); s = h.FontSize; h.FontSize = 12; j = h.Rotation; h.Rotation=45; k = h.Position; h.Position = [0.5 0.5 0]; b = ta1.Color; d=ta1.LineWidth; ta1.Color = 'red'; ta1.LineWidth= 4;
Semen: PotD and Spermidine Binding MATLAB Code
Before using MATLAB the original system of equations was non-dimensionalised to the system:
$$ \begin{eqnarray*} \frac{dP}{d\tau}&=&C-\kappa P S,\\ \frac{dS}{d\tau}&=&C-\kappa P S,\\ \frac{dC}{d\tau}&=&\kappa P S -C. \end{eqnarray*} $$The initial conditions were also non-dimensionalised to become:
$$ \begin{eqnarray*} P(0)&=&R_{0},\\ S(0)&=&1,\\ C(0)&=&0. \end{eqnarray*} $$To perform senstivity analysis that is described in the Semen section of the BioSpray models, MATLAB was used. Two files were written one to set the function, spermsen.m, and one to solve the function and plot the results in Figure 2, Figure 3 and Figure 4, run_spermsen.m. Both files are shown below, where the green is comments to aid in understanding of the scripts.
MATLAB code from spermsen.m file.
%Function to define the non-dimensionalised system of ODEs describing the %binding between Spermidine and PotD. u is a 3 dimensional vector %where; % u(1)= P (PotD Concentration), % u(2)= S (Spermidine Concentration) % u(3)= C (PotD and Spermidine Complex Concentration. ) % t is the time that the simulation is run over. % kappa is the binding parameter determoned by the binding rates and % initial concentration of Spermidine. function f = spermsen(t,u,kappa) % define the system of ODEs by setting the vector f = left hand side of the system. f = [u(3)-kappa.*u(1)*u(2); u(3)-kappa.*u(1)*u(2); kappa.*u(1)*u(2)-u(3)]; % The function is then ended. end % This function can then be called in the run_spermsen.m file by using % @spermsen command. %MATLAB code for run_spermsen.m file.
%%Section 1: File to perform sensitivity analysis for PotD and Spermidine binding model.This section runs the analysis on both the binding parameter kappa %and the ratio between PotD and Spermidine, u0. % a is the number of datapoints in the range for both kappa and u0, chosen to be 50. a=50; % kappa1 is the range of values for kappa that the system will be solved over, % where the max is twice the expected value. kappa1=linspace(0,(129.0875)*2,a); % u0 is the range of values for u0 that the system will be solved over % where the max is four times the expected value. u0=linspace(0,2,a); % T is the final time of the simulation in seconds. T=0.10; % Store is set as an empty matrix that will be filled by the values for the % complex concentration for each of the values of kappa and u0. Store = zeros(a,a); % The for loop solves the function spermsen.m using ode23 for the range of values of % kappa and u0, and stores the complex concentration values in the Store matrix. for i=1:a for j=1:a [t,u]=ode23(@spermsen,[0 T],[u0(a-i+1) 1 0],[],kappa1(j)); Store(i,j) = u(end,3); end end %imagesc plots the complex concentration against u0 and kappa. imagesc(Store) % label function adds x labels to both the x and y axis. xlabel('Increasing R_{0}','FontSize',15,'FontWeight','bold'); ylabel('Increasing \kappa','FontSize',15,'FontWeight','bold'); % the set function removes the values from the x and y axis. set(gca,'YTick',[],'XTick',[]); % the color bar function adds a user defined colourbar with the label, % 'Complex Concentration'. c=colorbar('Ticks',[0,0.99],'TickLabels',{'None','High'},'FontSize',15,'FontWeight','bold'); c.Label.String = 'Complex Concentration'; % The resulting plot is Figure 2. %%Section 2: This section is used to plot only u0 against complex concentration over time. a is the number of datapoints in the range for u0, chosen to be 50. a=50; % u0 is the range of values for u0 that the system will be solved over % where the max is four times the expected value. u0=linspace(0,2,a); % store is set as an empty matrix that will be filled by the values for the % time values, PotD concentration, Spermidine concentration and % complex concentration for each value of u0. store = zeros(a,4,a); %1st is time, 2:4 are u. % The for loop solves spermsen.m for the range of values of u0 and stores % tha values in the store matrix. for i=1:a [t,u]=ode23(@spermsen,[0 0.1],[u0(i) 1 0],[],129.0875); store(end-size(t,1)+1:end,1,i) = t; store(end-size(t,1)+1:end,2:4,i) = u; end % figure calls up an empty figure so that figure 2 is not overwritten. figure % u0 is then plotted against complex concentration, u(3). plot(u0,squeeze(store(end,4,:)),'LineWidth',2); % label function adds x labels to both the x and y axis. xlabel('Increasing R_{0}','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % This plot is Figure 3. %%Section 3: This section is used to plot only kappa against complex concentration over time. a is the number of datapoints in the range for kappa, chosen to be 50.a=50; % kappa1 is the range of values for kappa that the system will be solved over, % where the max is twice the expected value. kappa1=linspace(0,(129.0875)*2,a); % store is set as an empty matrix that will be filled by the values for the % time values, PotD concentration, Spermidine concentration and % complex concentration for each value of kappa. store = zeros(a,4,a); %1st is time, 2:4 are u % The for loop solves spermsen.m for the range of values of kappa and stores % tha values in the store matrix. for i=1:a [t,u]=ode23(@spermsen,[0 0.1],[1.4 1 0],[],kappa1(i)); store(end-size(t,1)+1:end,1,i) = t; store(end-size(t,1)+1:end,2:4,i) = u; end % figure calls up an empty figure so that figure 3 is not overwritten. figure % kappa is then plotted against complex concentration, u(3). plot(kappa1,squeeze(store(end,4,:)),'LineWidth',2); % axis tight ensures the axis are not too large for the datapoints. axis tight % label function adds x labels to both the x and y axis. xlabel('Increasing \kappa','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % This plot is Figure 4.
Saliva: Lactoferrin and Lactoferrin Binding Protein Binding MATLAB Code
Nasal Mucus: Oderant Binding Protein Folding MATLAB Code.
Before using MATLAB the original system of equations was non-dimensionalised to the system:
$$ \begin{eqnarray*} \frac{d\alpha}{d\tau}&=&OBP-\psi \alpha \beta,\\ \frac{d\beta}{d\tau}&=&OBP-\psi \alpha \beta,\\ \frac{dOBP}{d\tau}&=&\psi \alpha \beta -OBP. \end{eqnarray*} $$The initial conditions were also non-dimensionalised to become:
$$ \begin{eqnarray*} \alpha(0)&=&1,\\ \beta(0)&=&D_{0},\\ OBP(0)&=&0. \end{eqnarray*} $$To perform senstivity analysis that is described in the Nasal Mucus section of the BioSpray models, MATLAB was used. Two files were written one to set the function, obp.m, and one to solve the function and plot the results in Figure 5, Figure 6 and Figure 7, run_obp.m. Both files are shown below, where the green is comments to aid in understanding of the scripts.
-->MATLAB code from obp.m file
%Function to define the non-dimensionalised system of ODEs describing the %folding of the Oderant Binding Protein (OBP). u is a 3 dimensional vector %where; % u(1)= alpha (Alpha Helix Concentration), % u(2)= beta (Beta Barrel Concentration) % u(3)= OBP (Folded OBP Concentration. ) % t is the time that the simulation is run over. % K represents psi, the binding parameter determined by the binding rates and % initial concentration of alpha helices. function f = obp(t,u,K) % define the system of ODEs by setting the vector f = left hand side of the system. f = [u(3)-K*u(1)*u(2); u(3)-K*u(1)*u(2); -u(3)+K*u(1)*u(2)]; % The function is then ended. end % This function can then be called in the run_obp.m file by using % @obp command.
