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