Matlab code

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