From 2008.igem.org
function [omega_H_dim, H_dim, H_dim2, H_dim3]=model(omega_H, omega_AH, phi)
%This function takes in 3 dimensionless parameters (omega_H = holincontribution,
%omega_AH= antiholin contribution and phi = coupling contribution)
Hc=1000; % Critical Holin Concentration (1000 holin/cell)
H_dim = [];
omega_H_dim =[];
for j=.001:1:200
A = omega_H*j-omega_AH*.001+phi; % linear term representing the critical holin concentration
M = [1 -A -phi*omega_H*j]; % polynomial representing the critical holin concentration
B = max(roots(M)); % returns the larger of the two roots of M
H_dim = [H_dim B/Hc];
omega_H_dim = [omega_H_dim omega_H*j];
end
H_dim2 = [];
omega_H_dim2 =[];
for j=.001:1:200
A = omega_H*j-omega_AH+phi; % linear term representing the critical holin concentration
M = [1 -A -phi*omega_H*j]; % polynomial representing the critical holin concentration
B = max(roots(M)); % returns the larger of the two roots of M
H_dim2 = [H_dim2 B/Hc];
omega_H_dim2 = [omega_H_dim2 omega_H*j];
end
H_dim3 = [];
omega_H_dim3 =[];
for j=.001:1:200
A = omega_H*j-omega_AH*100+phi; % linear term representing the critical holin concentration
M = [1 -A -phi*omega_H*j]; % polynomial representing the critical holin concentration
B = max(roots(M)); % returns the larger of the two roots of M
H_dim3 = [H_dim3 B/Hc];
omega_H_dim3 = [omega_H_dim3 omega_H*j];
end
semilogx(50:1:5e3,ones(size(50:1:5e3)))
hold on
semilogx(omega_H_dim, H_dim, omega_H_dim2, H_dim2, omega_H_dim3, H_dim3)
axis([50, 5e3, 0.01, 2])
title(['Lambda Lysis Model with phi = ',int2str(phi)])
xlabel('omega H')
ylabel('[H]/[H]critical')
legend('0.001*omega AH','omega AH','100*omega AH')
hold off
"
