function main % This is MATLAB code for ODE yaws simulations % Authors Presley Kimball, Jacob Levenson, Amy Moore, Jan Rychtar, Dewey Taylor % Email questions to Jan Rychtar rychtarj@vcu.edu % Version from December 7, 2021 %----------- % Written on Matlab 2020a (version 9.8.0.1417392 (R2020a) Update 4') % with optimization toolbox % The program has several major parts, each starts with %% % % The main parts are as follows % % (Part 0) Technical preliminaries (Users do not have to touch) % % (Part 1) Setting up parameters (names and values) % This is a large part further divided into % 1a) List of all parameters and their base values % 1b) Technical part to set the parameters as empty variables to make them global. % 1c) Set auxiliary variables as functions of parameters % 1d) Actually set the parameter values to their base values % % (Part 2) Setting up compartments % % (Part 3) Setting up the equations for the system % % (Part 4) Actual code to run % % (Part 5) Codes of the functions that depend on the model % This is the actual core of the code. However, this part generally % does not have to be touched by a casual user that is interested % only in modifying a few parameter values % % This includes code for % - analytical solutions of the equilibria % - analytical solution of the basic reproduction number % - solving when R0=1 % - calibration % - uncertainty analysis % - validation % % (Part 6) Codes of functions that do not depend on model % This is a large part that essentially should not be modified by a % casual user % Functions in this part should be working for all sorts of models % - Resetting Parameters Functions % - Sensitivity analysis Functions % - Setting up system of equations to use for ODE solvers % - numerically solving the system of equations % - checking for consistency % % (Part 7) Some auxiliary functions for better plotting, easy saving etc % The important textual display is saved in the file yaws_code_output.txt %% (Part 0) Technical preliminaries % Close all figures close all; % Clear all variables clear all; % Set the randomg number generator for reproducibility rng(0) % Check if the folder for figures exists and create if not if not(isfolder('Figures')) mkdir('Figures') end if license('test','optimization_toolbox')~=1 disp_Diary('Optimization toolbox is not installed') disp_Diary('Calibration and uncertainty functions may not run properly') pause end % Set the file for textual output and write the date and time we run the code % After multiple run, the last output is at the bottom of this text file diary 'yaws_code_output.txt' disp_Diary('--------------- New run ----------------') disp_Diary(datetime) %% (Part 1) Setting up parameters (names and values) % Part 1a - list of all parameters, base values and ranges % Set the main parameters % Each line is for one parameter % The first column is the name that will be used in Matlab % The second column is the label that will be used for plotting % The third is the base value % The fourth is the minimal and the fifth is the maximal value (for % sensitivity analysis) Params = ... {'Lambda', '$\Lambda$', 27.2/12000, 0.001, 0.003; % Birth rate 'mu', '$\mu$', 1/(65*12), 0.0008, 0.0017; % Mortality rate 'beta', '$\beta$', 0.014273, 0.01, 0.02; % Transmission rate 'sigma', '$\sigma$', 30/21, 1/3, 10/3; % Incubation rate 'lambda1', '$\lambda_1$',1/3, 1/6, 1/2; % Healing rate of primary yaws 'lambda2', '$\lambda_2$',1/3, 1/6, 1/2; % Healing rate of secondary yaws 'rho1', '$\rho_1$', 1/1.5, 1/24, 1; % Primary latency period length 'rho2', '$\rho_2$', 1/25, 1/60, 1/12; % Secondary latency period length 'pL2Y2', '$p_{L_2Y_2}$', 0.9999, 0.9998, 1; % Probability of relapse from secondary latent to secondary yaws 'pY1Y2', '$p_{Y_1Y_2}$', 0.12, 0.09, 0.15; % Probability of primary yaws to secondary yaws 'pY2Y3', '$p_{Y_2Y_3}$', 0.0001, 0, 0.0002; % Probability of secondary yaws to tertiary yaws }; % Part 1b - Set the parameters as global variables % This enables them to be passed from a function to function without actually calling them % It may not be the cleanest way to do it from coding perspective, but % increases legibility of the code in some places % For functionality, we have to manual assign of all the variables here % once as empty Lambda = []; mu = []; epsilon = 0; % This one hardoced and is not going to change beta = []; sigma = []; lambda1 = []; lambda2 = []; rho1 = []; rho2 = []; pL2Y2 = []; pY1L1 = []; pY1Y2 = []; pY2L2 = []; pY2Y3 = []; % The ones below are for treatment rates cE = []; cY1 = []; cY2 = []; cY3 = []; cL1 = []; cL2 = []; tauM = 0; % Just initialize the treatment rate as 0 % Part 1c - Set the auxiliary variables and make them global parameters % (they are not parameters but are calculated from them easily and % they simplify the model) % Each line is for one aux variable % The first column is the name that will be used in Matlab % The second column is the label that will be used for plotting % The third is the formula AuxVars = {'tauE' , '$\tau_E$', 'cE .* tauM'; 'tau1' , '$\tau_{Y_1}$', 'cY1 .* tauM'; 'tau2' , '$\tau__{Y_2}$', 'cY2 .* tauM'; 'tau3' , '$\tau_{Y_3}$', 'cY3 .* tauM'; 'tauL1', '$\tau_{L_1}$', 'cL1 .* tauM'; 'tauL2', '$\tau_{L_2}$', 'cL2 .* tauM'; 'pY1L1', '$p_{Y_1L_1}$', '1-pY1Y2'; 'pY2L2', '$p_{Y_2L_2}$', '1-pY2Y3'; 'N', '$N$', 'Lambda./mu' }; % Set all auxiliary variables as global tauE = []; tau1 = []; tau2 = []; tau3 = []; tauL1 = []; tauL2 = []; N=[]; % Part 1d - Actually set the values of the parameters % Get the number of parameters, discard the number of columns [numParams,~] = size(Params); % % Get the number of auxiliary variables, discard the number of columns [numAuxVars,~] = size(AuxVars); % % Get the parameter names Names = Params(:,1); % Get the parameter labels Labels = Params(:,2); % Get the parameter base values BaseValues = Params(:,3); % Actually set the parameters ResetParameters % Tentatively set treatment as total community treatment TCT; %% (Part 2) Set up the compartments Comps = {'S', '$S$'; % Susceptible 'E', '$E$'; % Exposed 'Y1', '$Y_1$'; % Primary yaws 'Y2', '$Y_2$'; % Secondary yaws 'Y3', '$Y_3$'; % Tertiray yaws 'L1', '$L_1$'; % Primarily Latently infected 'L2', '$L_2$'; % Secondarily Latently infected }; % Set the compartments as empty global S = []; E = []; Y1 = []; Y2 = []; Y3 = []; L1 = []; L2 = []; dS = []; dE = []; dY1 = []; dY2 = []; dY3 = []; dL1 = []; dL2 = []; % Get the number of compartments, discard the number of columns [numComps,~] = size(Comps); % %% Part 3 set up the equations of the system function equations = getEquations % Set the equations % This returns the right hand sides of the equations % It is done via a function so that the equations are set when % needed. They can be used for symbolic solutions as well as for % numerical solutions or even for differential equations equations = { ... Lambda+tauE*E+(tau1 +(1-(pY1Y2+pY1L1))*lambda1)*Y1+(tau2+(1-(pY2Y3+pY2L2))*lambda2)*Y2+tau3*Y3+tauL1*L1+tauL2*L2-(epsilon+beta*(Y1/(S+E+Y1+Y2+Y3+L1+L2)+Y2/(S+E+Y1+Y2+Y3+L1+L2))+mu)*S; %dS/dt (epsilon+beta*(Y1/(S+E+Y1+Y2+Y3+L1+L2)+Y2/(S+E+Y1+Y2+Y3+L1+L2)))*S-(sigma+tauE+mu)*E; %dE/dt sigma*E-(lambda1+tau1+mu)*Y1; %dY1/dt pY1Y2*lambda1*Y1+rho1*L1+pL2Y2*rho2*L2-(lambda2+tau2+mu)*Y2; %dY2/dt pY2Y3*lambda2*Y2+(1-pL2Y2)*rho2*L2-(tau3+mu)*Y3; %dY3/dt pY1L1*lambda1*Y1-(rho1+tauL1+mu)*L1; %dL1/dt pY2L2*lambda2*Y2-(rho2+tauL2+mu)*L2; %dL2/dt }; end %% Part 4 Actual code to run starts here % Check for errors in the code/math % Run with an argument larger than 1 for more through check % It will generate random parameters specified number of times and % crosschecks the formulas and analytical solutions and R0 etc % checkConsistency(1) % Calibrate, i.e. find beta, and display the outcomes calibration('display') % Calculate the value of R0 disp(['R0 is ' num2str(R0)]) %pause % Produce the validation graphs Samples = 100; validation2 dependence_on_cL2 %Compute Critical values for tau for TCT and TTT %Create a plot of value of tau versus value of R0 plotR0onTau % Do the global stability globalStability(1000) % Do the uncertainty analysis % Set to 100 for good and 1000 for better % The function uncertainty saves the results, so make sure it runs again % with more samples by renaming the saved file Samples = 1000; % How many samples for sensitivity analysis % Actually set the parameters ResetParameters % Note that it generally takes a long time to run the uncertainty analysis % Thus, the procedure loads previously saved results and one should rename % or delete the file YawsUncertainty.mat if one want to see new pictures. uncertainty return %% (Part 5) Essential Functions that depend on the model function TTT(coverage) % Resets treatments probabilities for targeted treatments % This is the WORST CASE SCENARIO treating only the active yaws if nargin == 0 % if the function is called without an argument, % assume the coverage is 1 coverage = 1; end cE = coverage; cY1 = coverage; cY2 = coverage; cY3 = coverage; cL1 = coverage; %coverage; cL2 = 0.1; %0.01; CalculateAuxiliaryVariables; end function TCT(coverage) % Resets treatment probabilities for total community treatment if nargin == 0 % if the function is called without an argument, % assume the coverage is 1 coverage = 1; end cE = coverage; cY1 = coverage; cY2 = coverage; cY3 = coverage; cL1 = coverage; cL2 = coverage; CalculateAuxiliaryVariables; end function sols = getAnalyticalSolutions % Returns analytical solutions to the equilibrium % equations % This function is coded in full generality for nonzero epsilon % Define constants for simplification vE = sigma + tauE + mu; vY1 = lambda1 + tau1 + mu; vY2 = lambda2 + tau2 + mu; vY3 = mu + tau3; vL1 = rho1 + tauL1 + mu; vL2 = rho2 + tauL2 +mu; kY1 = sigma / vY1; kL1 = (pY1L1 * lambda1 * sigma) / (vY1 *vL1); kY2 = ((pY1Y2 * lambda1 * kY1 + rho1 * kL1)* vL2) / (vY2*vL2 - pL2Y2 * rho2 * pY2L2 * lambda2); kL2 = (pY2L2 * lambda2 * kY2) / vL2; kY3 = (pY2Y3 * lambda2 * kY2 + (1-pL2Y2)* rho2 *kL2) / vY3; kN = 1+kY1 +kY2 +kY3 + kL1 + kL2; w0 = -epsilon; w1= vE+ kN*epsilon -vE*R0; w2 = vE*kN*R0; % Solve the quadratic for E/N, pick only the positive solution EoverN = (-w1+sqrt((w1)^2-(4*w2*w0)))/(2*w2); %Positive solutions % Get E by mutliplying E/N by N ansE = EoverN * N; if ansE ~= 0 % If endemic equilibrium aS = vE * ansE/ (epsilon+beta*(kY1+kY2)*EoverN); else % Disease-free equilibrium aS = N; end % Calculate the rest of the compartments as constant times E aE = ansE; aY1 = kY1*ansE; aY2 = kY2*ansE; aY3 = kY3*ansE; aL1 = kL1*ansE; aL2 = kL2*ansE; % Return the analytical solutions sols = [aS, aE, aY1, aY2, aY3, aL1, aL2]; end function output = R0 % Returns the basic reproduction number for the system % Calculate the outgoing rates for the simplification vE = sigma + tauE + mu; vY1 = lambda1+tau1+mu; vY2 = lambda2+tau2+mu; vY3 = mu + tau3; vL1 = rho1+tauL1+mu; vL2 = rho2+tauL2 + mu; % Return R0 output = (beta.*sigma)./(vE.*vY1) + (beta.*(lambda1.*pY1L1.*rho1.*sigma.*vL2 + ... lambda1.*pY1Y2.*sigma.*vL1.*vL2))./(vE.*vL1.*(vL2.*vY1.*vY2 - lambda2.*pL2Y2.*pY2L2.*rho2.*vY1)); end function output = calibration(input) % Calibrates the model by finding the values of beta, pY1L1 and % pY2L2 to fit data from Lihir Island % "input" is an optional parameters. If any is given, it will % display the results of the calibration. % output is the error after the calibration function output = getPrevalences(params) % Calculates prevalences of S, active and latent in the % population, % params give beta % output is the prevalences % Set paramters as given beta = params(1); CalculateAuxiliaryVariables; % Get analytical solutions sol = getAnalyticalSolutions; % Output the prevalences output = [sol(1)+sol(2), sol(3)+sol(4)+sol(5), sol(6)+sol(7)]/sum(sol); end % Do the optimization only if toolbox is installed if license('test','optimization_toolbox')==1 % Data from observations Mitja et al (2018) realData = [0.797, 0.024, 0.179]; % Set the function to minimize (it is the norm distance of our estimated % prevalences from the real data) fun = @(params) norm((realData - getPrevalences(params))./realData); % Set the lower bound constraint (rates and probabilities are non-negative) lb = [0]; % Set the upper bound constraint (beta probabilities have to add to no more than 1) ub = [1]; % Set initial guess init = [0]; % There are no linear or nonlinear constraints, so set those arguments to []. A = []; b = []; Aeq = []; beq = []; nonlcon = []; % Suppress the display of minimization results options = optimoptions('fmincon','Display','off'); % Do the actual minimization subject to the constraints specified above x = fmincon(fun,init,A,b,Aeq,beq,lb,ub, nonlcon, options); else disp('Optimization toolbox not installed') disp('Using hardcoded values') x = beta; end % Unpack the results beta = x; % Calculate the error of the approximation output = fun(beta); if nargin>0 disp(['Calibration complete. Optimal values are ' ... ' beta = ' num2str(beta) ]); disp(['Real data are ' num2str(realData)]) disp(['Our estimated values are ' num2str(getPrevalences(x))]) disp(['Error is ' num2str(output)]) end end function results = solve_for_tau % Numerically solves where R0=1 and outputs the solution % Store what the original value of tauM is tauOriginal = tauM; % Check whether there is a solution or not tauM = 0; % sets treatment to 0 CalculateAuxiliaryVariables % this sets all treatment rates appropriate values if R0<1 results = 0; % R0<1 even without treatment, so the solution is 0 else % Create a symbolic parameter for tau TAU = sym('TAU'); % Set tau as a symbolic variable tauM = TAU; CalculateAuxiliaryVariables % Set the equation as R0=1. eqn = (R0 == 1); % Numerically solve the equation % There will be 5 solutions, only one is positive and the maximal % is the right solution tauSolve = vpasolve(eqn, tauM); results = tauSolve(tauSolve>0); end % Restore the original value of tau tauM = tauOriginal; CalculateAuxiliaryVariables; end function results = solve_for_tau_TTT % Numerically solves where R0=1 and outputs the solution under TTT TTT; results = solve_for_tau; end function dependence_on_cL2 tauM = 0; % Get init conditions based on calibrated model (i.e. close to % current true data) ResetParameters; CalculateAuxiliaryVariables; InitConditions = getAnalyticalSolutions; all_cL2_values = linspace(0.01,1,100); for aux1 = 1:length(all_cL2_values) % Start with a year of TCT treatment tspan = [0, 12]; % 1 year tauM = 1/6; TCT; % Numerically solve the ODE [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, InitConditions); % Get prevalence of all infections as a function of time I_TCT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time) % Get the times needed for elimination decrease = 1000; % Now we will simulate TTT after two rounds of TCT % See where we stopped after a year (two rounds of TCT) Init_afterTCT = y(end,:); % Reset the time to run again for 40 years tspan = [0, 480]; % Set the treatment rate to reflect the rate from the % validation tauM = 1/6; % Set the value of cL2 cL2 = all_cL2_values(aux1); CalculateAuxiliaryVariables; % Numerically solve the ODE [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, Init_afterTCT); % Get prevalence of all infections as a function of time I_TTT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time) % Get the truly initial conditions I_TTT(1) = 1-InitConditions(1)/(sum(InitConditions)); % Get how long it will take to decrease the prevalence given number of times if I_TTT(end)R0>1 and % fitting error is small while error(aux1)>0.005 %(R0Un(aux1)<=1)||(R0Un(aux1)>3)|| % Assign parameter values at random (but within the bounds) for aux = 1:numParams eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';']) end % Set tau = 0 for no treatment tauM = 0; % Calibrate to find beta and qY1 and qY2, store fittting % error for next round of "while cycle" error(aux1) = calibration; % Calculate the auxiliary variables CalculateAuxiliaryVariables % Evaluate R0 R0Un(aux1) = R0; end % Now we have random parameters such that 112, 1),:); % Reset the time to run again for 40 years tspan = [0, 480]; % Set the treatment rate to reflect the rate from the % validation tauM = 1/6; %tauM = 1/3;1/3 TTT; % Numerically solve the ODE [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, Init_afterTCT); % Get prevalence of all infections as a function of time I_TTT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time) % Get the truly initial conditions I_TTT(1) = 1-InitConditions(1)/(sum(InitConditions)); if aux1 == 1 figure plot(t,I_TTT) xlabel('Time (in months)') ylabel('Proportion of infected') FormatFigure; end % Get how long it will take to decrease the prevalence given number of times if I_TTT(end) 20 save(savedFileName,'R0Un','paramsUn', 'SolsUn', 'tElimTCT', 'tElimTTT', 'R0treatment', 'tEst', 'error' ) end end % Process and display the results of the uncertainty analysis disp(['Time to eliminate by TCT with base values ' num2str(tElimTCT(1)/12) ' years']) disp(['Probability to eliminate in 5 years ' num2str(mean((tElimTCT<60)))]) disp(['Probability to eliminate in 10 years ' num2str(mean((tElimTCT<120)))]) disp(['Time to eliminate by TTT with base values ' num2str(tElimTTT(1)/12) ' years']) disp(['Probability to eliminate in 10 years ' num2str(mean((tElimTTT<120)))]) disp(['Probability to eliminate in 15 years ' num2str(mean((tElimTTT<15*12)))]) disp(['Probability to eliminate in 20 years ' num2str(mean((tElimTTT<240)))]) disp(['Probability to eliminate in 25 years ' num2str(mean((tElimTTT<300)))]) % Investigate what is going on with sigma (nothing significant % found) % for aux = 1:numParams % figure % sc = scatter(paramsUn(:,4), paramsUn(:,aux), ... % 'MarkerFaceColor', [0.5 0.5 0.5], ... % 'MarkerEdgeColor', [0.5 0.5 0.5]); % alpha(sc, 1/Samples) % sc.SizeData = 5; % xlabel(Labels{4}) % ylabel(Labels{aux}) % % Set limits on y axis based on parameter limits % ylim([Params{aux,4} Params{aux,5}]) % FormatFigure % SaveFigure(['Correlation_sigma_' Names{aux}]) % % % end % % See how well we can estimate time to elimination from R0 % % (NOT AS WELL AS WE HOPED FOR) % figure % sc = scatter(tElimTTT/12, tEst, ... % 'MarkerFaceColor', [0.5 0.5 0.5], ... % 'MarkerEdgeColor', [0.5 0.5 0.5]); % alpha(sc, 1/Samples) % sc.SizeData = 5; % xlabel('Time to elimination (in years)') % ylabel('$-(\log(R_0))^{-1}$') % FormatFigure % SaveFigure('Correlation between times') figure %edges = [10:1:31]; % histogram(tElimTTT/12, edges, 'Normalization', 'probability'); histogram(tElimTTT/12,'Normalization', 'probability'); xlabel('Time of TTT application (in years)') ylabel('Probability of elimination') FormatFigure SaveFigure(['tElimTTT_hist_Marino']) % Plot the bars with PRCC indices prcc = PRCClocal(paramsUn, tElimTTT', 'tElimTTT') figure histogram(tElimTCT/12, 'Normalization', 'probability'); xlabel('Time of TCT application (in years)') ylabel('Probability of elimination') FormatFigure SaveFigure(['tElimTCT_hist_Marino']) % Plot the bars with PRCC indices prcc = PRCClocal(paramsUn, tElimTCT', 'tElimTCT') pause figure histogram(R0treatment, 'Normalization', 'probability'); xlabel('$R_0$ with TTT treatment') ylabel('Probability') FormatFigure SaveFigure(['R0Treatment_hist_Marino']) % Plot the bars with PRCC indices prcc = PRCClocal(paramsUn, R0treatment', 'R0treatment') figure histogram(R0Un, 'Normalization', 'probability'); xlabel('$R_0$ prior treatment') ylabel('Probability') FormatFigure SaveFigure(['R0_hist_Marino']) % Plot the bars with PRCC indices prcc = PRCClocal(paramsUn, R0Un', 'R0') % figure % histogram(error, 'Normalization', 'probability'); % xlabel('Fit error') % ylabel('Probability') % FormatFigure % SaveFigure(['error_hist_Marino']) % % % Plot the bars with PRCC indices % prcc = PRCClocal(paramsUn, error', 'error') % % % Plot correlation between sigma and compartments % for i1 = 1:numComps % figure % sc = scatter(paramsUn(:,4), SolsUn(:,i1), ... % 'MarkerFaceColor', [0.5 0.5 0.5], ... % 'MarkerEdgeColor', [0.5 0.5 0.5]); % alpha(sc, 1/Samples) % xlabel(Labels{4}) % ylabel(Comps{i1,2}) % % FormatFigure % SaveFigure(['Correlations_sigma_' Comps{i1,1}]) % % cor = corrcoef(paramsUn(:,4),SolsUn(:,i1)) % disp(['Correlation between sigma and ' Comps{i1,1} ' is ' num2str(cor(1,2))]) % end % Plot correlations between compartments for i1 = 1:numComps for i2 = 1:numComps figure if i1==i2 histogram(SolsUn(:,i1), 'Normalization', 'probability'); xlabel(Comps{i1,2}) ylabel('Probability') else sc = scatter(SolsUn(:,i1), SolsUn(:,i2), ... 'MarkerFaceColor', [0.5 0.5 0.5], ... 'MarkerEdgeColor', [0.5 0.5 0.5]); alpha(sc, 1/Samples) xlabel(Comps{i1,2}) ylabel(Comps{i2,2}) end FormatFigure SaveFigure(['Correlations' num2str(i1) num2str(i2)]) end end % Plot parameter distribution for i1 = 1:numParams figure histogram(paramsUn(:,i1), 'Normalization', 'probability'); xlabel(Labels{i1}) ylabel('Probability') FormatFigure SaveFigure([Names{i1} '_hist']) end end function plotR0onTau % Plots how R0 depends on tau and graphically solves the equation % R0=1 % Store what tau was tauOrig = tauM; tauMax = 0.005; % Reset tau for plotting tauM = linspace(0,tauMax, 101); % Calculate R0 during TTT TTT; % Sets probabilities c to TTT setting R0TTT = R0; % Stores the R0 values tau0TTT = solve_for_tau; % Solves for tau0 during TTT dispTTT = ['Tau for TTT is ', num2str(sym2poly(tau0TTT)), ' which corresponds to a round of TTT every ', num2str(1/(12*sym2poly(tau0TTT))), ' years.']; disp(dispTTT) % Calculate R0 during TCT TCT; % Sets probabilities c to TTT setting R0TCT = R0; % Stores the R0 values tau0TCT = solve_for_tau; % Solves for tau0 during TCT dispTCT = ['Tau for TCT is ', num2str(sym2poly(tau0TCT)), ' which corresponds to a round of TCT every ', num2str(1/(12*sym2poly(tau0TCT))), ' years.']; disp(dispTCT) % Plot R0 on tau figure pTTT = plot(tauM, R0TTT, 'k--', 'DisplayName', 'TTT'); hold on pTCT = plot(tauM, R0TCT, 'k-', 'DisplayName', 'TCT'); % Plot horizontal line where R0 is 1 plot(tauM, ones(length(tauM),1), ':', 'color', [0.5 0.5 0.5]) % Plot vertical lines where tau = tau 0 plot([tau0TTT, tau0TTT], [0,1], ':', 'color', [0.5 0.5 0.5]) plot([tau0TCT, tau0TCT], [0,1], ':', 'color', [0.5 0.5 0.5]) % Set x and y labels xlabel('Treatment rate $\tau$') ylabel('Basic reproduction number') % Show the legend lgd = legend([pTTT pTCT], 'Location', 'NorthEast', 'Interpreter', 'LaTeX'); % Specify x axis labels xticks([0 sym2poly(tau0TCT) 0.001 sym2poly(tau0TTT) 0.002 0.003 0.004 0.005 ]) xticklabels({'','$\tau_{0,TCT}$', '', '$\tau_{0,TTT}$', '', '0.003', '0.004', '0.005'}) yticks([0 0.5 1 1.253 1.5]) xlim([min(tauM) max(tauM)]) % Format and save the figure FormatFigure SaveFigure('R0ontau') % Reset parameters back to what they were ResetParameters CalculateAuxiliaryVariables end function check_TTT % This checks whether analytical solutions for TTT are correct tauM = 0; % sets treatment to 0 CalculateAuxiliaryVariables % this sets all treatment rates appropriate values %Defining constants for simplification r1 = (lambda1*(rho2+mu))/(rho1+mu); r2= pY1L1*rho1+pY1Y2*(rho1+mu); r3= (rho2+mu)*(lambda2+mu)-pY2L2*lambda2*pL2Y2*rho2; r4 = (sigma + mu)*(lambda1 + mu); r5 = r3*r4*R0; t2= (rho2+mu)*(sigma+mu)*cY1*cY2; t1 = (rho2+mu)*cY2*r4+(sigma+mu)*cY1*r3-beta*sigma*(rho2+mu)*cY2; t0=r3*r4-r5; %Calculate Tau critical display((-t1+sqrt(t1^2-4*t2*t0))/(2*t2)); %display((-t1-sqrt(t1^2-4*t2*t0))/(2*t2)); % Reset the parameters ResetParameters CalculateAuxiliaryVariables end function validation % Fits real data to the curves and produces the figure ResetParameters %calibration CalculateAuxiliaryVariables % Set the data % Number of cases of active yaws (Table 1 from % https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5920722/ ) month = [0 6 12 18 24 30 36 42]'; cases = [323/13490 44/13166 34/13204 38/15977 24/11792 52/14935 53/14765 63/13601]'; latent = [18.9 14.7 7.3 10.8 5.4 7.9 4.7 6.6]'/100; % Plot the results % It will be done in two figures at the same time, one for active % and one for latent % Start the active infections activeFigure = figure; axActive = gca; % plot active infections real data plot(axActive, month,cases, '.', 'MarkerSize', 15); hold on % Start the latent infections latentFigure = figure; axLatent = gca; plot(axLatent, month-3,latent, '.', 'MarkerSize', 15); hold on % Plot our model prediction % Plots the solutions of the ODE system (specified by the equations listed in the function getEquations) % Specify the time of interest, first x months of total community % treatment tspan = [0 6]; % Specify initial conditions as endemic equilibrium InitConditions = getAnalyticalSolutions/N; % One round of TCT, i.e. the rate is 1/time tauM = 1/max(tspan); TCT; % Numerically solve the ODE. % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely) % The input ODEmodel is a function that has to be defined separately [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions); % Unpack the solutions (give it human names) % sols is an array: the first component is time (stored in t), the second component are the compartments. for comp = 1:numComps % Evaluate the value of the compartment eval([Comps{comp} ' = sols(:,comp) ;']); % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution % but for all compartments automatically and in one line end % Plot the active infections in the proper axis (figure) plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k'); % Plot the latent infections in the proper axis (figure) plot(axLatent, t,(L1+L2)./(S+E+L1+L2), 'k') % Specify the time of interest, up to 24 months of TTT tspan = [max(tspan) 24]; % Start where we ended InitConditions = sols(end,:); % 3 rounds of TTT plus additional non-strategic treatment between % rounds tauM = 1/6; %tauM = 1/3; TTT; % Numerically solve the ODE. % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely) % The input ODEmodel is a function that has to be defined separately [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions); % Unpack the solutions (give it human names) % sols is an array: the first component is time (stored in t), the second component are the compartments. for comp = 1:numComps % Evaluate the value of the compartment eval([Comps{comp} ' = sols(:,comp) ;']); % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution % but for all compartments automatically and in one line end plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k--'); plot(axLatent, t, (L1+L2)./(S+E+L1+L2), 'k--'); % Specify the time of interest tspan = [24 42]; % Start where we ended InitConditions = sols(end,:); %Set treatment rate - Non-strategic treatment tauM = 1/24; TTT; %Numerically solve the ODE. %ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely) %The input ODEmodel is a function that has to be defined separately [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions); %Unpack the solutions (give it human names) %sols is an array: the first component is time (stored in t), the second component are the compartments. for comp = 1:numComps %Evaluate the value of the compartment eval([Comps{comp} ' = sols(:,comp) ;']); %The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution %but for all compartments automatically and in one line end plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k:'); plot(axLatent, t,(L1+L2)./(S+E+L1+L2), 'k:'); %Set axis ticks, labels, and bounds xticks(axLatent, month) xlabel(axLatent, 'Month') ylabel(axLatent, 'Prevalence of latent yaws') xlim([-6,42]) ylim([0, 0.2]) % Format the figure FormatFigure(latentFigure) % Save the figure SaveFigure('Validation_Latent_Yaws_TCT_TTT_lowerTTT') close(latentFigure) %Set axis ticks, labels, and bounds xticks(axActive, month) xlabel(axActive, 'Month') ylabel(axActive, 'Prevalence of active yaws') xlim([0,42]) ylim([0, 0.03]) % Format the figure FormatFigure(activeFigure) % Save the figure SaveFigure('Validation_Active_Yaws_TCT_TTT_lowerTTT') % Display output for the caption of the figure disp('Validation output to the caption') tauM = 0; disp(['$R_0$ Predicted by the model is ' num2str(R0)]); end function validation2 % Fits real data to the curves and produces the figure with some % sort of uncertainty ResetParameters %calibration CalculateAuxiliaryVariables % Set the data % Number of cases of active yaws (Table 1 from % https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5920722/ ) month = [0 6 12 18 24 30 36 42]'; cases = [323/13490 44/13166 34/13204 38/15977 24/11792 52/14935 53/14765 63/13601]'; latent = [18.9 14.7 7.3 10.8 5.4 7.9 4.7 6.6]'/100; myColor = [0 0 0]; % Sets the original color as black % Plot the results % It will be done in two figures at the same time, one for active % and one for latent % Start the active infections activeFigure = figure; axActive = gca; % plot active infections real data plot(axActive, month,cases, '.', 'MarkerSize', 15); hold on % Start the latent infections latentFigure = figure; axLatent = gca; plot(axLatent, month-3,latent, '.', 'MarkerSize', 15); hold on tauM = 0; % Get init conditions based on calibrated model (i.e. close to % current true data) ResetParameters; CalculateAuxiliaryVariables; InitConditions = getAnalyticalSolutions; % For the number of experiments, generate parameters randomly and % see what the outcomes are for aux1 = 1:Samples+1 disp(['Running simulation ' num2str(aux1) ' of ' num2str(Samples) ]) if aux1 == 1 % Use the base values for the first run % We will just set R0 as it is and then the values will not % have to be generated again R0Un(aux1) = R0; error(aux1) = 0; else % Tentatively set R0 to 0 so that we can discard parameters % that produce "bad" R0 and too large error R0Un(aux1) = 0; error(aux1) = 1; end % Keep assigning parameter values randomly until 3>R0>1 and % fitting error is small while error(aux1)>0.9 %(R0Un(aux1)<=1)||(R0Un(aux1)>3)|| % Assign parameter values at random (but within the bounds) for aux = 1:numParams eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';']) end % Set tau = 0 for no treatment tauM = 0; % Calibrate to find beta and qY1 and qY2, store fittting % error for next round of "while cycle" error(aux1) = calibration; % Calculate the auxiliary variables CalculateAuxiliaryVariables % Evaluate R0 R0Un(aux1) = R0; end % Now we have random parameters such that 11 corresponds to endemic equilibrium (by % comparing it to numerical solutions % ADD MORE AS NEEDED disp('Checking for consistency. It may take some time') % Initialize the number of times R0>1 Endemic = 0; % Initialize the number of time endemic is unstable Unstable = 0; for dummy_variable= 1:num_of_tests disp(['Test ' num2str(dummy_variable) ' of ' num2str(num_of_tests) '.' ]) % Assign parameter values at random (but within the specified range) for aux = 1:numParams % Go through all the parameters % Evaluate the command that assigns the parameter with a given name the proper base value %eval([Names{aux} '=' num2str(BaseValues{aux}*2*rand) ';']) eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';']) end % Calculate the auxiliary variables CalculateAuxiliaryVariables % Calculate Analytical Solutions ansols = getAnalyticalSolutions % Check whether analytical solutions are correct tol = 10^(-10); % tolerance level for correctness of analytical solutions if max(abs(cell2mat(CheckSolution(ansols))))>tol % if the max error for analytical solutions is bigger then % tolerance, display the warning disp('Analytical solutions may be wrong.') displayParameters end % Check whether R0<1 corresponds to no endemic equilibrium % Get R0 Reff = R0; % Update the number of times R0>1 if Reff>1 Endemic = Endemic + 1; end % Get numerical solutions numsols = numSolve; if (((Reff<1) && (any(numsols.E>0))) && (epsilon ==0)) % if R0<1 but there is an endemic equilibrium (positive exposed) disp('R0<1 but endemic equilibrium exists') displayParameters end if ((Reff>1) && (all(numsols.E<=0))) % if R0>1 but there is no endemic equilibrium (no positive exposed) disp('R0>1 but endemic equilibrium does not exist') displayParameters end n_of_numsols = []; eval(['n_of_numsols = size(numsols.' Comps{1} ',1);']); % Get the index where numerical and analytical should agree solIdx = 0; % Initialize the index for aux = 1:n_of_numsols % For every numerical solution if abs(numsols.E(aux) - ansols(2))1 % only if the endemic equilibrium is possible disp('Checking global stability of endemic equilibria') % Set initial conditions InitConditions = [1 1 0 0 0 0 0]; % Set thenumber of months to run the ODEs tspan = [0, 12000]; % Numerically solve the ODE [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, InitConditions); % Get the difference betweennumerical and analytical % solution difference = max(abs(ansols - y(end,:))); if difference >10^(-3) disp('Equilibrium may not be stable') Unstable = Unstable +1; displayParameters end end end disp(['Out of ' num2str(Endemic) ' endemic equilibria, ' num2str(Unstable) ' were unstable']) disp('Consistency check complete. If no warning was displayed, no inconsistency was found') % Reset parameters back to base values ResetParameters CalculateAuxiliaryVariables end function prcc=PRCClocal(values, Y, FileName); % PRCC procedure % input: Y, response function that we investigate how it depends % on the parameters % output - the prcc coefficients % Based on Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008). % A methodology for performing global uncertainty and sensitivity analysis in systems % biology. Journal of Theoretical Biology, 254(1):178–196 % Part of the code for PRCC taken from % http://malthus.micro.med.umich.edu/lab/usanalysis.html % Initiate the LHS matrix that will carry all the parameters LHSmatrix= values; for i=1:numParams % Loop for the parameters % Remove the ith parameter LHStemp=LHSmatrix; LHStemp(:,i)=[]; Ztemp=LHStemp; LHStemp=[]; % Get the correlation coefficient [rho_aux,~]=partialcorr([LHSmatrix(:,i),Y],Ztemp,'type','Spearman'); % Record the correlation coefficient prcc(1,i)=rho_aux(1,2); end % Sort the correlation coefficients for better plotting [sortedPRCC,I] = sort(prcc, 'descend'); % % Plot PRCC coefficients as a horizontal bar plot % figure; % h=barh(sortedPRCC); % % Set the y limits by the number of parameters % ylim([0.5 numParams+0.5]); % % Ticks at integers and labeled by the parameter names/labels % ax = gca; % ax.YTick = 1:numParams; % ax.YTickLabel = Labels(I); % ylabel('Parameters'); % % Set the x axis between -1 and 1 % xlim([-1 1]); % xlabel('PRCC'); % % Make the color gray % h.FaceColor = [0.5 0.5 0.5]; % % Format the figure % FormatFigure; % SaveFigure([FileName '_bar_Marino']) % % PLot only the important coefficients figure % Determine which ones are important importantPRCC = sortedPRCC(abs(sortedPRCC)>0.01); h=barh(importantPRCC); ylim([0.5 length(importantPRCC)+0.5]); % Ticks at integers and labeled by the parameter names/labels ax = gca; ax.YTick = 1:length(importantPRCC); ax.YTickLabel = Labels(I(abs(sortedPRCC)>0.01)); ylabel('Parameters'); % Set the x axis between -1 and 1 xlim([-1 1]); xlabel('PRCC'); % Make the color gray h.FaceColor = [0.5 0.5 0.5]; FormatFigure; SaveFigure([FileName '_bar_Marino_important']) end function out = displayParameters % Displays parameter values for cross-checking disp('Parameter values are') for aux = 1:numParams % Go through all the parameters and display their values eval([Names{aux}]) end pause end function globalStability(numOfRepetition) % Checks whether the endemic equilibrium is globally stable % Defines a Lyapunov function and checks that the derivative is % negative by randomly generating a bunch of compartments and % pluggin it in ResetParameters; % Store what the original value of tauM is tauOriginal = tauM; % Check whether there is a solution or not tauM = 0; % sets treatment to 0 TCT(0); % Get analytical solutions for endemic equilibrium asols = getAnalyticalSolutions; % Unpack the solutions aS = asols(1); aE = asols(2); aY1 = asols(3); aY2 = asols(4); aY3 = asols(5); aL1 = asols(6); aL2 = asols(7); for aux = 1:numOfRepetition % Randomly assign values to all compartments for comp = 1:numComps % for all compartments % Assign a value to the compartment eval([Comps{comp,1} '=' num2str(rand(1)*10) ';']) end % Evaluate the derivatices ders = cell2mat(getEquations); % Unpack the derivatives for comp = 1:numComps % for all compartments % Assign a value to the compartment eval(['d' Comps{comp,1} '=' num2str(ders(comp)) ';']) end % Evaluate the derivative to Lyapunov function dLyap(aux) = (1-aS/S)*dS + (1-aE/E)*dE +(1-aY1/Y1)*dY1 +... (1-aY2/Y2)*dY2 +(1-aY3/Y3)*dY3 +(1-aL1/L1)*dL1 + ... (1-aL2/L2)*dL2; end % Check that no value is positive if max(dLyap)>0 disp('SOMETHING IS WRING WITH GLOBAL STABILITY') else disp('Endemic equilibrium seems globally stable') end % Reset tauM to whatever it was tauM = tauOriginal; end %% (Part 7) Auxiliary Functions for easy plotting, etc function FormatFigure(h) % Formats figures so that all figures are uniformly formatted if nargin == 0 % if the function is called without an argument, assume we want to format current figure h = gcf; end % Get axes handle haxes=get(h, 'CurrentAxes'); % Get figure number fnbr = get(h,'Number'); [fig_x,fig_y] = getPosition(fnbr); % Set dimensions of the figure set(h,'Units','Inches','Position',[fig_x, fig_y, 3, 2.5]); % Set default font size set(haxes,'FontSize',10) % Set interpreter to LaTeX set(0,'defaultTextInterpreter','latex') % Set the interpreter for latex set(haxes, 'TickLabelInterpreter', 'latex'); % Find all line objects hline= findobj(haxes,'Type','line'); % Set line properties (relatively fat lines) set(hline, 'LineWidth' , 1); % Set axes properties set(haxes, ... 'Box' , 'on' , ... %puts a box around the axis 'TickDir' , 'in' , ... 'TickLength' , [.02 .01] , ... 'XMinorTick' , 'off' , ... 'YMinorTick' , 'off' , ... 'LineWidth' , 1); function [xpos,ypos] = getPosition(n) % gets the x y position of the figure to plot n = mod(n-1,15)+1; ypos = floor((n-1)/5)*5; xpos = (mod(n-1,5))*3; end end function SaveFigure(filename) % Saves a current graphics file to file named filename.pdf % Get the handle h=gcf; set(h,'Units','Inches'); pos = get(h,'Position'); % Cut all the white margins set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]) % Actually save print(h,['Figures/' filename, '.pdf'],'-dpdf','-r0') % This outputs on the screen what figure was saved into what file disp_Diary(['Saved a figure ' num2str(get(gcf,'Number')) ' into a file ' ... filename '.pdf']) end function disp_Diary(string) % Displays text on the screen and also logs it in the diary diary on disp(string) diary off end end