Team:Dundee/Modeling/Appendix1
To see the MATLAB code for the Fingerprint Ageing section of the project or the Chromate Biosensor section of the project use the following buttons.
Appendix 3: Chromate Biosensor Appendix 2: Fingerprint Ageing
Blood: MATLAB Code
The figure shown in the haptoglobin and haemoglobin binding model were created using MATLAB. 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*} $$Two files were written to perform sensitivity analysis; 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)=\alpha (free haemoglobin concentration), % u(3)= [\alpha \cdot Hp] (haptoglobin and haemoglobin initial complex concentration), % u(4)= [\alpha \cdot Hp \cdot \beta] (final haemoglobin/haptoglobin complex concentration). % Parameter t is the time that the simulation is run over. % Parameter lambda represents the parameter in the system determined by binding rates and initial concentration of haemoglobin. % Parameter 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 = right 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 the @haptosen command. % END OF FILE MATLAB code from run_haptosen.m file % File to perform sensitivity analysis for haptoglobin and haemoglobin binding model. % Variable a is the number of values chosen for each of the ranges of parameters, 100 different values were chosen.. a=100; % Define lambda1 as the range of values for the parameter lambda, chosen with the maximum value as twice the expected value. lambda1=linspace(0,((100/317)*0.3878494524)*2,a); % Define gamma1 as the range of values for the parameter gamma, chosen with the maximum value as twice the expected value. gamma1=linspace(0,(83/317)*2,a); % T is the final time of the simulation, in seconds. T=30; % Define Store as an empty matrix for the final concentrations of the complex with the varied parameter values. Store = zeros(a,a); % Command 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); end end % Command 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'); % Command set removes the x and y axis values. set(gca,'YTick',[],'XTick',[]); % The following lines add a colour bar with defined label 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; % END OF FILE↑ View description of model
Semen: MATLAB Code
The figures shown in the spermidine and potD binding model were created using MATLAB. 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*} $$Two files were written to perform sensitivity analysis; 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. % Parameter 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). % Parameter t is the time that the simulation is run over. % Parameter kappa is the binding parameter determined 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 = right 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 the @spermsen command. % END OF FILE 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 (R0 in the described model). % Define a as the number of datapoints in the range for both kappa and u0, chosen to be 50. a=50; % Define kappa1 as the range of values for kappa that the system will be solved over, where the maximum value is twice the expected value. kappa1=linspace(0,(129.0875)*2,a); % Define u0 as the range of values for u0 that the system will be solved over, where the maximum value 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 % Command imagesc plots the complex concentration against u0 and kappa. imagesc(Store) % Command label 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. % Define a as the number of datapoints in the range for u0, chosen to be 50. a=50; % Define u0 as the range of values for u0 that the system will be solved over, where the maximum value is four times the expected value. u0=linspace(0,2,a); % Define store 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 the 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 % Command figure calls up an empty figure so that figure 2 is not overwritten. figure % Then u0 can be plotted against complex concentration, u(3). plot(u0,squeeze(store(end,4,:)),'LineWidth',2); % Command label 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. % Define a as the number of datapoints in the range for kappa, chosen to be 50. a=50; % Define kappa1 as the range of values for kappa that the system will be solved over, where the maximum value is twice the expected value. kappa1=linspace(0,(129.0875)*2,a); % Define store 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 the 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 % Command figure calls up an empty figure so that figure 3 is not overwritten. figure % Then kappa can be plotted against complex concentration, u(3). plot(kappa1,squeeze(store(end,4,:)),'LineWidth',2); % Command axis tight ensures the axis are not too large for the datapoints. axis tight % Command label 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. % END OF FILE↑ View description of model
Saliva: MATLAB Code
The figures shown in the lactoferrin and LBP binding model were created using MATLAB. Before using MATLAB the original system of equations was non-dimensionalised to the system:
$$ \begin{eqnarray*} \frac{dLf}{d\tau}&=&[Lf \cdot LBP]-\theta Lf LBP,\\ \frac{dLBP}{d\tau}&=&[Lf \cdot LBP]-\theta Lf LBP,\\ \frac{d[Lf \cdot LBP]}{d\tau}&=&\theta Lf LBP -[Lf \cdot LBP]. \end{eqnarray*} $$The initial conditions were also non-dimensionalised to become:
$$ \begin{eqnarray*} Lf(0)&=&1,\\ LBP(0)&=&v_{0},\\ [Lf \cdot LBP](0)&=&0. \end{eqnarray*} $$Two files were written to perform sensitivity analysis; one to set the function (lactofun.m) and one to solve the function and plot the results in Figure 5, Figure 6 and Figure 7 (run_lactofun.m). Both files are shown below, where the green is comments to aid in understanding of the scripts.
MATLAB code for lactofun.m file. % Function to define the non-dimensionalised system of ODEs describing the binding between lactoferrin and LBP. % Parameter u is a 3 dimensional vector where; % u(1) represents lactoferrin, % u(2) represents the lactoferrin binding protein, % u(3) represents the final complex. % Parameter t is the time that the simulation is run over. % Parameter theta is the binding parameter determined by the binding rates and the initial concentration of lactoferrin. function f = lactofun(t,u,theta) % Define the system of ODEs by setting the vector f = right hand side of the system. f=[u(3)-theta.*u(1)*u(2); u(3)-theta.*u(1)*u(2); -u(3)+theta.*u(1)*u(2)]; % The function is then ended. end % This function can then be called in the run_lactofun.m file by using % @lactofun command. %END OF FILE. MATLAB code for run_lactofun.m file. %% Section 1: File to perform sensitivity analysis for lactoferrin and lactoferrin binding protein binding model. %This section runs the analysis on both the binding parameter theta and the ratio between lactoferrin binding protein and lactoferrin, v0. % Variable a is the number of datapoints in the range for both theta and v0, chosen to be 50. a=50; % Define theta1 as the range of values for theta that the system will be solved over, where the maximum value is twice the expected value. theta1=linspace(0,623.417*2,a); % Define v0 as the range of values for v0 that the system will be solved over where the maximum value is four times the expected value. v0=linspace(0,4,a); % T is the final time of the simulation in seconds. T=10; % Set Store as an empty matrix that will be filled by the values for the complex concentration for each of the values of theta and v0. Store = zeros(a,a); % The for loop solves the function lactofun.m using ode23 for the range of values of theta and v0, and stores the complex concentration values in the Store matrix. for i=1:a for j=1:a [t,u]=ode23(@lacto,[0 T],[1 v0(a-i+1) 0],[],theta1(j)); Store(i,j) = u(end,3); end end % Call up a new figure. figure(1) % Command imagesc plots the complex concentration against v0 and theta. imagesc(Store) % The label function adds x labels to both the x and y axis. xlabel('Increasing v_{0}','FontSize',15,'FontWeight','bold'); ylabel('Increasing \theta','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 5. % % %% Section 2: This section is used to plot only v0 against complex concentration over time. % Define store1 as an empty matrix that will be filled by the values for the time values, lactoferrin concentration, lactoferrin binding protein concentration and complex concentration for each value of v0. store1 = zeros(5000,4,a); % 1st is time, 2:4 are u. % The for loop solves lactofun.m for the range of values of v0 and stores the values in the store1 matrix for i=1:a [t,u]=ode23(@lactofun,[0 T],[1 v0(i) 0],[],5); store1(end-size(t,1)+1:end,1,i) = t; store1(end-size(t,1)+1:end,2:4,i) = u; end % Command figure calls up an empty figure so that Figure 5 is not overwritten. figure (2) % Then v0 can be plotted against complex concentration, u(3). plot(v0,squeeze(store1(end,4,:)),'LineWidth',2); % The label function adds x labels to both the x and y axis. xlabel('Increasing v_{0}','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % This plot is Figure 6. % % %% Section 3: This section is used to plot only theta against complex concentration over time. % Define store2 as an empty matrix that will be filled by the values for the time values, lactoferrin concentration, lactoferrin binding protein concentration and complex concentration for each value of theta. store2 = zeros(5000,4,a); % 1st is time, 2:4 are u % The for loop solves lactofun.m for the range of values of theta and stores the values in the store2 matrix. for i=1:a [t,u]=ode23(@lactofun,[0 T],[1 1 0],[],theta1(i)); store2(end-size(t,1)+1:end,1,i) = t; store2(end-size(t,1)+1:end,2:4,i) = u; end % Command figure calls up an empty figure so that Figure 6 is not overwritten. figure(3) % Then theta can be plotted against complex concentration, u(3). plot(theta1,squeeze(store2(end,4,:)),'LineWidth',2); % Command axis defines the axis we wish to consider. axis([0,1200,0,1]) % The label function adds x labels to both the x and y axis. xlabel('Increasing \theta','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % This plot is Figure 7. % END OF FILE.↑ View description of model
Nasal Mucus: MATLAB Code.
The figures shown in the odorant binding protein folding model were created using MATLAB. 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*} $$Two files were written to perform sensitivity analysis; one to set the function (obp.m) and one to solve the function and plot the results in Figure 8, Figure 9 and Figure 10 (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 odorant 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). % Parameter t is the time that the simulation is run over. % Parameter 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 = right 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 the @obp command. % END OF FILE MATLAB code for run_obp.m file. % File to perform sensitivity analysis for OBP folding model. %% Section 1: This section runs sensitivity analysis for both the parameter psi (K) and the ratio between beta barrels and alpha helices, v0 (D0 in the described model). % Define a as the number of datapoints in the range for both K and v0, chosen to be 50. a=50; % Define K as the range of values for K that the system will be solved over, where the maximum value is twice the expected value. K=linspace(0,(5291.005292)*2,a); % Define u0 as the range of values for v0 that the system will be solved over, where the maximum value is four times the expected value. v0=linspace(0,4,a); % T is the final time of the simulation in seconds. T=0.001; % Define Store as an empty matrix that will be filled by the values for the folded OBP concentration for each of the values of K and v0. Store = zeros(a,a); % Call up first empty figure to plot in. figure(1) % The for loop solves the function obp.m using ode23 for the range of values of K and v0, and stores the complex concentration values in the Store matrix. for i=1:a for j=1:a [t,u]=ode23(@obp,[0 T],[1 v0(a-i+1) 0],[],K(j)); Store(i,j) = u(end,3); end end % The command imagesc plots the complex concentration against v0 and K. imagesc(Store) % The label function adds x labels to both the x and y axis. xlabel('Increasing D_{0}','FontSize',15,'FontWeight','bold'); ylabel('Increasing \psi','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 8. % % %% Section 2: This section is used to plot only v0 against complex concentration over time. a is the number of datapoints in the range for v0, chosen to be 50. a=50; % Define v0 as the range of values for v0 that the system will be solved over, where the maximum value is four times the expected value v0=linspace(0,4,a); % define store as an empty matrix that will be filled by the values for the time values, alpha helix concentration, beta-barrel concentration and folded OBP concentration for each value of v0. store = zeros(a,4,a); % 1st is time, 2:4 are u. % The for loop solves obp.m for the range of values of v0 and stores the values in the store matrix. for i=1:a [t,u]=ode23(@spermsen,[0 0.001],[1 v0(i) 0],[],5291.005292); store(end-size(t,1)+1:end,1,i) = t; store(end-size(t,1)+1:end,2:4,i) = u; end % The command figure calls up an empty figure so that Figure 8 is not overwritten. figure(2) % Then v0 can be plotted against complex concentration, u(3). plot(v0,squeeze(store(end,4,:)),'LineWidth',2); % The label function adds x labels to both the x and y axis. xlabel('Ratio of \beta to \alpha, D_{0}','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % The annotation command allows us to add lines with specified style, width and colour for annotation of the graph, these can then be manually manipulated by using the edit plot tool. t2 = annotation('line', [0.42 0.42], [0.15 0.91]); b = t2.Color; d=t2.LineWidth; g=t2.LineStyle; t2.Color = 'red'; t2.LineWidth= 2; t2.LineStyle= ':'; % This plot is Figure 9. % % %% Section 3: This section is used to plot only K against folded OBP over time. % Define a as the number of datapoints in the range for K, chosen to be 50. a=50; % Define K as the range of values for K that the system will be solved over, where the maximum value is twice the expected value. K=linspace(0,(5291.005292)*2,a); % Define as an empty matrix that will be filled by the values for the time values, alpha helix concentration, beta-barrel concentration and folded OBP concentration for each value of K. store = zeros(a,4,a); % 1st is time, 2:4 are u for i=1:a [t,u]=ode23(@spermsen,[0 0.001],[1.5 1 0],[],K(i)); store(end-size(t,1)+1:end,1,i) = t; store(end-size(t,1)+1:end,2:4,i) = u; end % The for loop solves obp.m for the range of values of K and stores the values in the store matrix. % The command figure calls up an empty figure so that Figure 9 is not overwritten. figure(3) % Then K can be plotted against complex concentration, u(3). plot(K,squeeze(store(end,4,:)),'LineWidth',2); % The command axis tight ensures the axis are not too large for the datapoints. axis tight % The label function adds x labels to both the x and y axis. xlabel('Increasing K','FontSize',15,'FontWeight','bold'); ylabel('Concentration of Complex Formed','FontSize',15,'FontWeight','bold'); % This plot is Figure 10. % END OF FILE↑ View description of model
