% % Example8_1.m % clear all global coef Y0 wt % % Set parameters and initial conditions % % Note: We split the coefficients between params, which is passed to ode45 % and coef which is global. The former comprises the set employed in % Chapters 12 and 13 for Bayesian inference and uncertainty propagation. % lam1 = 1e+4; d1 = .01; epsilon = 0; k1 = 8.0e-7; lam2 = 31.98; d2 = 0.01; f = 0.34; k2 = 1e-4; delta = 0.7; m1 = 1.0e-5; m2 = 1.0e-5; NT = 100; c = 13; rho1 = 1; rho2 = 1; lamE = 1; bE = 0.3; Kb = 100; dE = 0.25; Kd = 500; deltaE = 0.1; coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE,dE,deltaE,m1,m2,Kd]; parameters = [bE; delta; d1; k2; lam1; Kb]; n = 41; p = 6; Y0 = [0.9e+6; 4000; 1e-1; 1e-1; 1; 12]; % % Compute and plot the Effector cell population on the interval [0,200] to illustrate % the dynamics for the nominal parameter values % t_data = 0:.1:200; ode_options = odeset('RelTol',1e-6); [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters); E = Y(:,6); figure(1) plot(t_data,E,'-k','linewidth',3) hold on plot(t_data(501),E(501),'ob','linewidth',5) hold off set(gca,'Fontsize',[22]); xlabel('Time (Days)') ylabel('Immune Effector Response E') % % Compute the Effector cell population at time T = 50 for the values in lamE_vec % and rho1 % per = 0.1; N_pert = 20; lamE_low = lamE-lamE*per; lamE_high = lamE+lamE*per; lamE_step = (lamE_high-lamE_low)/N_pert; lamE_vec = lamE_low:lamE_step:lamE_high; rho1_low = rho1-rho1*per; rho1_high = rho1+rho1*per; rho1_step = (rho1_high-rho1_low)/N_pert; rho1_vec = rho1_low:rho1_step:rho1_high; tvals = 0:.01:50; N = length(tvals); for j=1:N_pert+1 lamE = lamE_vec(j); coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE,dE,deltaE,m1,m2,Kd]; parameters = [bE; delta; d1; k2; lam1; Kb]; [t,Y] = ode15s(@hiv_rhs,tvals,Y0,ode_options,parameters); V_lamE(j) = Y(N,5); E_lamE(j) = Y(N,6); end lamE = 1; for j=1:N_pert+1 rho1 = rho1_vec(j); coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE,dE,deltaE,m1,m2,Kd]; parameters = [bE; delta; d1; k2; lam1; Kb]; [t,Y] = ode15s(@hiv_rhs,tvals,Y0,ode_options,parameters); V_rho1(j) = Y(N,5); E_rho1(j) = Y(N,6); end figure(2) plot(rho1_vec,E_rho1,'-k',lamE_vec,E_lamE,'--b','linewidth',3); axis([0.9 1.1 34 42]) set(gca,'Fontsize',[23]); xlabel('Parameter Value') ylabel('Immune Effector Response E') legend('\rho_1','\lambda_E','location','SouthEast') figure(3) plot(rho1_vec,V_rho1,'-k',lamE_vec,V_lamE,'--b','linewidth',3); axis([0.9 1.1 1.525e+5 1.55e+5]) set(gca,'Fontsize',[23]); xlabel('Parameter Value') ylabel('Free Virus V') legend('\rho_1','\lambda_E','location','SouthEast') % % End hiv_model %