MATLAB code from squalene.m file
% File to solve the system of ODEs describing lanosterol synthase and squalene epoxide binding.
% The vector u has 4 dimensions;
% u(1) = Lanosterol synthase (LS) concentration, {non-dimensionalized}
% u(2) = Squalene epoxide (SE) concentration, {"}
% u(3) = Protosterol cation (PC) concentration (1st Intermediate), {"}
% u(4) = Lanosterol (La) concentration. {"}
% t is the time that the simulation is run over.
% lambda is the 1st binding parameter ((k1/k2)SE(0)).
% gamma is the 2nd binding parameter (k3/k2).
function f = squalene(t,u,lambda,gamma);
% Then set f as the RHS of the system of ODEs.
f = [u(3)-lambda.*u(1)*u(2);
    u(3)-lambda.*u(1)*u(2);
    lambda.*u(1)*u(2)-u(3)-gamma.*u(3);
    gamma.*u(3)];
end
% Then end the function.
% This function can then be solved using ode23 in the file run_squalene.m  by using the @squalene command.
% END OF FILE.


MATLAB code for run_squalene.m file.
% File to solve the non-dimensionalised system of ODEs representing the binding of squalene epoxide (SE) and lanosterol synthase (LS) to form lanosterol (La).
%
% Set the final time of the simulation, set at 30 seconds.
T = 30;
% Set initial concentration of ratio (LS/SE) as the expected value.
u0 = 2.561;
% Parameter a is how many values of each parameter are used in the sensitivity analysis.
a = 100;

% Set lambda1 as the range of values for parameter lambda, where the maximum value is twice the expected.
lambda1 = linspace(0,(0.64061499*2),a);
% Set gamma1 as the range of values for parameter gamma, where the maximum value is twice the expected.
gamma1 = linspace(0,(1.000223733*2),a);
% Set Store as an empty matrix to store the values of lanosterol concentrations for each range of values of lambda1 and gamma1.
Store = zeros(a,a);
% Set a for loop that will solve the system of ODEs from the squalene.m file using the ode23 solver, for each value of lambda1, gamma1.

% The resulting values for the lanosterol concentration are stored in the Store matrix.
for i=1:a
    for j=1:a
[t,u]=ode23(@squalene,[0 T],[u0 1 0 0],[],lambda1(a+1-i),gamma1(j));
 Store(i,j) = u(end,4);

    end
end
% Call up a blank figure, so that no other figures are overwritten.
figure(1)
% Plot a colourplot of Store, that is plot the lanosterol concentration against lambda1 and gamma1.

imagesc(Store)
% Set the colormap as default colour choice, as this is easiest to see.
colormap default
% Remove the values on the x and y axis by using the set command.
set(gca,'YTick',[],'XTick',[]);
% Add a colorbar with user defined axis labels (Tick labels) and label,
% 'Complex Concentration'.
c=colorbar('Ticks',[0,0.99],'TickLabels',{'None','High'},'FontSize',15,'FontWeight','bold');
c.Label.String = 'Complex Concentration';
% Add labels to the x and y axis for clarity on the plot.
xlabel('\fontsize{12} \bf Increasing \lambda');
ylabel('\fontsize{12} \bf Increasing \gamma');
%
%
%% Section 2: This section is used to plot only lambda against lanosterol concentration over time.
% Define store1 as an empty matrix to hold values of LS, SE, PC (1st Intermediate) and La, for each value of lambda1.
store1 = zeros(5000,5,a); % 1st is time, 2:5 are u.
% Set a for loop that will solve the system of ODEs from the squalene.m file using the ode23 solver, for each value of lambda1, with the expected value of gamma1 used.
% The resulting values for LS, SE, PC and the lanosterol concentration are stored in the store1 matrix.
for i=1:a

    [t,u]=ode23(@squalene,[0 T],[2.561 1 0 0],[],lambda1(i),1.000223733);
    store1(end-size(t,1)+1:end,1,i) = t;
    store1(end-size(t,1)+1:end,2:5,i) = u;

end
% Call up a blank figure, so that no other figures are overwritten.
figure(2)
% Plot the values of lambda1 against the corresponding values of lanosterol concentration over time.

 plot(lambda1,squeeze(store1(end,5,:)),'LineWidth',2);
 % Add labels to the x and y axis for clarity on the plot.
 xlabel('\bf \fontsize{12}Increasing \lambda');
 ylabel('\bf \fontsize{12} Lanosterol Concentration');
 % Set the axis to tight so that the axis is not too large for the data.
 axis tight
 %
 %
%% Section 3: This section is used to plot only gamma against lanosterol concentration over time.
% Define store2 as an empty matrix to hold values of LS, SE, PC and La, for each value of gamma1.
store2 = zeros(a,5,a); % 1st is time, 2:4 are u.
% Set a for loop that will solve the system of ODEs from the squalene.m file using the ode23 solver, for each value of gamma1, with the expected value of lambda1 used.
% The resulting values for LS, SE, PC and the lanosterol concentration are stored in the store2 matrix.
for i=1:a

    [t,u]=ode23(@squalene,[0 T],[2.561 1 0 0],[],0.64061499,gamma1(i));
    store2(end-size(t,1)+1:end,1,i) = t;
    store2(end-size(t,1)+1:end,2:5,i) = u;

end
% Call up a blank figure, so that no other figures are overwritten.
figure(3)
% Plot the values of gamma1 against the corresponding values of lanosterol concentration over time.
 plot(gamma1,squeeze(store2(end,5,:)),'LineWidth',2);
% Add labels to the x and y axis for clarity on the plot.
  xlabel('\bf Increasing \gamma');
 ylabel('\bf Lanosterol Concentration');
% Set the axis to tight so that the axis is not too large for the data.
 axis tight
 %
 %
%% Section 4: This section is used to plot only v0, the ratio of (LS(0)/SE(0)) against lanosterol concentration over time.
% Define store3 as an empty matrix to hold values of LS,SE,PC and La, for each value of v0.
store3= zeros(a,5,a);
% Set v0 as the range of values for v0 with the maximum as twice the expected value.
v0=linspace(0,(2.561*2),a);
% Set a for loop that will solve the system of ODEs from the squalene.m file using the ode23 solver, for each value of v0, with the expected value of lambda and gamma used.
% The resulting values for LS, SE, PC and the lanosterol concentration are stored in the store3 matrix.
for i=1:a
    [t,u]= ode23(@squalene,[0 T],[v0(i) 1 0 0],[],0.64061499,1.00223733);
    store3(end-size(t,1)+1:end,1,i) = t;
    store3(end-size(t,1)+1:end,2:5,i) = u;
end
% Call up a new figure, so that no figures are overwritten.
figure(4)
% Plot the values of v0 against the corresponding values of lanosterol concentration over time.
plot(v0,squeeze(store3(end,5,:)),'LineWidth',2);
% Set the axis to tight so that the axis is not too large for the data.
axis tight;
% Add labels to the x and y axis for clarity on the plot.
xlabel('\bf \fontsize{12} Increasing v_{0}');
ylabel('\bf \fontsize{12} Concentration of Lanosterol formed');
%
% END OF FILE.