CODEWITHSUNDEEP
matlabode
materialb0y
materialb0y

Reputation: 21

How to make bifurcation plot for the genetic toggle switch model?

I wanted to make a bifurcation plot to demonstrate the bistability region for the following system of equations: Toggle switch equations

I have already determined the nullclines for this system as well as their intersection points (via MATLAB). I read from here that I'll have to plug in the said points in the Jacobian, but I honestly have no idea how to proceed from that. For the bifurcation plot, I wanted to plot alpha1 vs alpha2.

Would appreciate any additional hints for the MATLAB implementation!

% This code block sets up the differential equations
% and the equations for the nullclines

alpha1 = 5;     % Effective synthesis rate of repressor 1
alpha2 = 5;     % Effective synthesis rate of repressor 2
beta = 3;       % Repression cooperativity of promoter 2
gamma = 3;      % Repression cooperativity of promoter 1

maxY = 5.5;                     %  Specifies max y-axis value

u = linspace(0,maxY,100);
v = linspace(0,maxY,100);

nullU = alpha1./(1+v.^beta);    % Setting du/dt = 0
nullV = alpha2./(1+u.^gamma);   % Setting dv/dt = 0

% I'm also exploring the time-course of u and v with respect to time
% so I'm using this function as well

function hillfx = toggle(t,y,alpha1,alpha2,beta,gamma)

    du = -y(1) + alpha1/(1+(y(2)^beta));
    dv = -y(2) + alpha2/(1+(y(1)^gamma));
    
    hillfx = [du; dv];
end

Upvotes: 2

Views: 153

Answers (1)

BoilingT
BoilingT

Reputation: 38

You could try varying 1 of the parameters. A good way to do this is to make a new variable for the range of alpha1 for example (this can be done with the linspace() function). Next create a matrix to store the equilibrium points from the resulting differential equations and set the initial conditions for the ode45 function (y0) to [0,0]. You will also need to make variables for beta and gamma accordingly.

Loop i over the length of alpha1Range and first save alpha1Range(i) to a variable which could be called a1. Use the ode45 function which solves differential equations with a1, alpha2, beta,gamma and y0. You might need to use the toggle function aswell. Save the output(time steps and solution values) of the ODE's for the current value of alpha1 to variables and then store the equilibrium points for this data in the i position of equilibrium points matrix.

You can now plot the bifurcation diagram by calling the plot function with alpha1Range and all the elements in the first column of equilibium points.

Upvotes: 1

Related Questions