% MIT Software licence % Copyright 2023 Rychtar and Taylor % % Permission is hereby granted, free of charge, to any person obtaining % a copy of this software and associated documentation files (the “Software”), % to deal in the Software without restriction, including without limitation % the rights to use, copy, modify, merge, publish, distribute, sublicense, % and/or sell copies of the Software, and to permit persons to whom % the Software is furnished to do so, subject to the following conditions: % % The above copyright notice and this permission notice shall be included % in all copies or substantial portions of the Software. % % THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS % OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, % FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL % THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER % LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, % OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS % IN THE SOFTWARE. % % function Schistosomiasis % This is MATLAB code for Math380 project on Schistosomiasis % modified by Santiago Q. Lopez, Samiya F. Majid and Rida A. Syed, Jan Rychtar and Dewey Taylor % based on the template by Jan Rychtar and Dewey Taylor % Email questions to Jan Rychtar rychtarj@vcu.edu % Version from August 16, 2023 % Written on Matlab 2022a (version 9.12.0.2009381) % The program has several major parts, each starts with %% % % (Part 0) Technical preliminaries (no need to touch) % % (Part 1) Setting up the compartments % (UPDATE ACCORDING TO YOUR MODEL) % % (Part 2) Setting up the parameters (names and values) % (UPDATE ACCORDING TO YOUR MODEL) % This is a large part further divided into % 2a) List of all parameters and their base values % Also contains list of controls % 2b) List of auxiliary variables as functions of parameters % 2c) Actually set the parameter values to their base values % % % (Part 3) Setting up the equations % (UPDATE ACCORDING TO YOUR MODEL) % % (Part 4) Actual code to run % (UPDATE ACCORDING TO WHAT YOU WANT TO DO WITH YOUR MODEL) % % (Part 5) Codes of the functions that depend on the model % (UPDATE ACCORDING TO YOUR CALCULATIONS) % 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 % - checking for consistency % - May contain calibration, validation and uncertainty analysis % % % (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 % - finding roots of functions % % (Part 7) Some auxiliary functions for better plotting, easy saving etc % The important textual output is saved in the file SIRV_code_output.txt %% (Part 0) Technical preliminaries % Close all figures close all; % Clear all variables clearvars; % Set the random number generator for reproducibility rng(0) % Check if the folder for figures exists and create if not if not(isfolder('Figures')) mkdir('Figures') 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 'Schisto.txt' disp_Diary('--------------- New run ----------------') disp_Diary(datetime) %% (Part 1) Set up the compartments Comps = {'S1', '$S_1$'; % Susceptible Humans 'I1', '$I_1$'; % Infected Humans 'V1', '$V_1$'; % Vaccinated humans 'S2', '$S_2$'; %Suseptable Snails 'I2', '$I_2$';% Infected Snails 'M', '$M$'; %Mericadia 'P', '$P$'; %Cericae }; % Set the compartments as empty ``global'' variables S1 = []; I1 = []; V1 = []; M = []; S2 = []; I2=[]; P=[]; % Get the number of compartments, discard the number of columns [numComps,~] = size(Comps); % %% (Part 2) Setting up parameters (names and values) % Part 2a - 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 = ... {'lambda1', '$\Lambda_1$', 0.031, 0.02, 0.04; % rate at which susceptible humans are recruited 'mu1', '$\mu_1$', (1/20), 1/25, 1/15; % Natural mortality rate for humans 'mu3', '$\mu_3$', 365/0.25, 1100, 1750; % Natural mortality rate for miracidia 'mu2', '$\mu_2$', 1.85, 1.5, 2.4; %natural death rate of snails 'mu4', '$\mu_4$', 830, 500, 1100; %natural death rate of the cericaie 'gamma2', '$\gamma_2$', 1.55*10^5, 0.9*10^5, 2.5*10^5; %rate of cericae released by infected snails 'gamma1', '$\gamma_1$', 1.1*10^5, 10^5, 2*10^5; % Miracidia hatch rate 'lambda2', '$\Lambda_2$', 10, 5, 15; %recruitment of susceptible snails 'delta1', '$\delta_1$', 1/10000, 0, 0.001; % Disease-related deaths of infected humans 'delta2', '$\delta_2$', 0.25, 0.2, 0.3; %Disease-related deaths for snail population 'beta1', '$\beta_1$', 0.0013, 10^(-3), 1.5* 10^(-3); % Probability of transmission on contact 'beta2', '$\beta_2$', 12.71, 10, 15; %probability transmission coefficient for snail 'alpha1', '$\alpha_1$', 0.0315, 0.01, 0.05; % scaling factor 'M0', '$M_0$', 3.5*10^3,3*10^3,5*10^3 %initial population of mircadia 'epsilon', ' $\epsilon$', 1.689, 1, 2; 'nu', '$\nu$', 0, 0, 0.1; % Vaccination rate 'omega', '$\omega$', 1/6.5, 1/8, 1/5; % vaccine waning 'CV', '$c$', 0.02, 0.001, 0.1; 'theta', '$\theta$', 0, 0, 0.1; %elimination rate of snails 'tau', '$\tau$', 0.0, 0.0, 0.1; %elimination rate of cercariae 'eta', '$\eta$', 0.0, 0, 0.2; %treatment rate of infected population }; % The controls follow the same structure as the parameters % They are used typically in vector borne diseases when people spray % against mosquitos, etc. Controls = { 'dummy', '$d$', 0, 0, 1; %elimination rate of snails }; % By default, set the controls to 0 (even if their values are not zero) % This will results in simulations as if no controls are used % We can change this setting anytime later setcontrols = 0; % Part 2b - 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 lambda1 = []; mu1 = []; mu3 = []; beta1 =[]; delta1= []; gamma1= []; M0= []; eta=[]; tau=[]; lambda2=[]; mu2=[]; mu4=[]; beta2=[]; delta2=[]; gamma2=[]; alpha1=[]; theta=[]; nu = []; omega = []; epsilon = []; dummy =[]; CD = 1; CV = []; % unit cost of vaccination % Part 2c - Set the auxiliary variables and make them global parameters % (they are not parameters but are calculated from them easily and % they simplify the model) % A typical use: if some parameters are given as times (duration of % infectious period) but the model needs the rate (1/duration), % we can address it here. % % Not so clean, but still working, usage is to have a total population size % calculated here from the parameters (works only if there is no disease % related mortality) % % 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 = { ... 'gamma', '$gamma$', '(lambda1*gamma1)/M0'; 'd1', '$d_1$', 'mu1+delta1+eta'; 'd2','$d_2$', 'mu2+theta'; 'd3', '$d_3$', 'mu2+delta2+theta'; 'd4', '$d_4$', 'mu4+tau'; 'delta', '$delta$','gamma2*lambda2'; 'alpha2', '$\alpha_2$', 'epsilon/M0'; }; % Set all auxiliary variables as global gamma = []; d1=[]; d2=[]; d3=[]; d4=[]; delta=[]; alpha2 = []; % Part 2d - Actually set the values of the parameters % Get the number of parameters, discard the number of columns [numParams,~] = size(Params); % [numControls,~] = size(Controls); % % Get the number of auxiliary variables, discard the number of columns [numAuxVars,~] = size(AuxVars); % % Get the parameter names Names = Params(:,1); % Get the parameter base values BaseValues = Params(:,3); % Actually set the parameters ResetParameters CalculateAuxiliaryVariables %% Part 3 set up the equations of the system function equations = getEquations % This returns the right hand sides of the equations associated % with the model % It is done via a function so that the equations are set whenever % needed. They can be used for symbolic solutions as well as for % numerical solutions or even for differential equations % The equations have to be in the same order as the compartments % (if the compartments are SIRV, the equations have to correspond % to S first, I second, R third and V last equations = { ... 1 - ((beta1.*P.*S1)/(1+(alpha1.*P)))-(mu1.*S1)+(eta.*I1)-(nu.*S1)+(omega.*V1); % dS1/dt ((beta1.*P.*S1)/(1+(alpha1.*P)))- d1.*I1; %dI1/dt nu.*S1-(omega+mu1).*V1; %dV1/dt gamma.*I1 - mu3.*M; %dM/dt 1-((beta2.*M.*S2)/(1+(alpha2.*M^2)))-(d2.*S2); %dS2/dt ((beta2.*M.*S2)/(1+(alpha2.*M^2)))-(d3.*I2);%dI2/dt (delta.*I2)-(d4.*P); %dP/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) % Set the parameters to their base values ResetParameters CalculateAuxiliaryVariables %justTry %EstimateParameterValues %return % Calculate the value of R0 and herd immunity R0 = getR0; disp(['R0 is ' num2str(R0)]) disp(['Herd imunity achieved at nu= ' num2str(getnuHI)]) ResetParameters CalculateAuxiliaryVariables disp(['Nash equilibrium achieved at p= ' num2str(getpNE)]) %Return the prevalence when everyone uses the nash equilibrium getPRne getPR(getnuHI) % Plot the figure that shows the calculation of herd imunity graphically plotCriticalValues(0,0.5,'nu', 'Vaccination rate, $\nu$') FormatFigure SaveFigure('Fig2') ResetParameters CalculateAuxiliaryVariables % Plot figure 3, the dependence of NE on the cost of protection plotFig3 ResetParameters CalculateAuxiliaryVariables % Plot Fig 4 - the dependence on eta plotDependenceOnEta ResetParameters CalculateAuxiliaryVariables % Plot sensitivities - figure 5,6,7 % Sensitivity of R, HI and NE (without controls) figure('Units','Inches','Position',[0, 0, 10, 5.5]); t=tiledlayout(2,3) t.Padding = 'compact'; t.TileSpacing = 'compact'; mytile = 1; LHS_PRCC(@() min(getR0,5), 'Reproduction number, $R_0$', 'R0', Params(1:18,:), 1000, mytile, mytile+3) mytile = mytile+1; LHS_PRCC(@() min(getnuHI,2), '$\nu_{HI}$', 'nuHI', Params(1:18,:), 1000, mytile, mytile+3) mytile = mytile+1; LHS_PRCC(@() min(getpNE,2), '$\nu_{NE}$', 'nuNE', Params(1:18,:), 1000, mytile, mytile+3); mytile = mytile+1; SaveFigure('Fig5') function output = getnuHImax1 % returns nuHI but truncated by 1 % this is to avoid very rare cases when nuHI goes beyond 1 % used for histogram plotting only output = min(getnuHI,1); end; function output = getpNEmax1 % returns nuHI but truncated by 1 % this is to avoid very rare cases when nuHI goes beyond 1 % used for histogram plotting only output = min(getpNE,1); end; % Sensitivity of the prevalence without controls figure('Units','Inches','Position',[0, 0, 6.6, 2.5]); t=tiledlayout(1,2) t.Padding = 'compact'; t.TileSpacing = 'compact'; mytile = 1; LHS_PRCC(@getPRne, 'Prevalence', 'PR', Params(1:18,:), 1000, mytile, mytile+1) SaveFigure('Fig6') % Sensitivity of the prevalence with controls figure('Units','Inches','Position',[0, 0, 6.6, 2.5]); t=tiledlayout(1,2) t.Padding = 'compact'; t.TileSpacing = 'compact'; mytile = 1; LHS_PRCC(@getPRne, 'Prevalence', 'PR', Params, 1000, mytile, mytile+1) SaveFigure('Fig7') return function plotFig3 % plots figure 3, the dependence of NE on the cost of protection % Get minimal value of K0 % (find which parameter is K0 by strcmpi) % minimum is the 4th column of the params variable % need to convert the cell to number at the end K0min = cell2mat( Params(strcmpi(Names,'CV'),4)); % Do the same to find the max K0max = cell2mat( Params(strcmpi(Names,'CV'),5)); % Set all possible values of K0 K = linspace(K0min,K0max,1001); for aux = 1:length(K) % for all values of K0 % Set the value of K0 CV = K(aux); % Calculate NE cNE(aux) = getpNE; % Calculate R0 at NE R(aux) = getR0vaccination(cNE(aux)); PR(aux) = getPR(cNE(aux)); end figure('Units','Inches','Position',[0, 0, 10, 2.5]); t=tiledlayout(1,3) t.Padding = 'compact'; t.TileSpacing = 'compact'; nexttile(1) plot(K,cNE, 'k-') xlabel('Relative cost of vaccination, $c$') ylabel('Voluntary vaccination, $\nu_{NE}$') FormatSubFigure min(K(cNE==0)) %SaveFigure('fig2') nexttile(2) plot(K,R, 'k-') %xticks(0:0.2:1) xlabel('Relative cost of vaccination, $c$') ylabel('Reproduction number, $R_0$') FormatSubFigure nexttile(3) plot(K,PR, 'k-') hold on plot([0 0.02], [0.01 0.01], 'k:') Kcrit = max(K(PR<0.01)) plot([Kcrit Kcrit], [0 0.01], 'k:') xlim([0 0.02]) %xticks(0:0.2:1) xlabel('Relative cost of vaccination, $c$') ylabel('Prevalence of infections') FormatSubFigure SaveFigure('Fig3') end function plotDependenceOnEta % Plots dependance of eta % Set all possible values of K0 alleta = linspace(0,0.1,101); for aux = 1:length(alleta) % for all values of K0 % Set the value of K0 eta = alleta(aux); CalculateAuxiliaryVariables % Calculate NE cNE(aux) = getpNE; nuHI(aux) = getnuHI; % Calculate R0 at NE PR(aux) = getPR(cNE(aux)); end figure('Units','Inches','Position',[0, 0, 6.5, 2.5]); t=tiledlayout(1,2) t.Padding = 'compact'; t.TileSpacing = 'compact'; nexttile(1) plot(alleta,cNE, 'k:') hold on plot(alleta, nuHI, 'k-') xlabel('Rate of MDA, $\eta$') ylabel('Vaccination rate') FormatSubFigure min(alleta(cNE==0)) min(alleta(nuHI==0)) %SaveFigure('fig2') nexttile(2) plot(alleta,PR, 'k-') hold on plot([0 0.1], [0.01 0.01], 'k:') Kcrit = min(alleta(PR<0.01)) plot([Kcrit Kcrit], [0 0.01], 'k:') xlim([0 0.1]) %xticks(0:0.2:1) xlabel('Rate of MDA, $\eta$') ylabel('Prevalence of infections') FormatSubFigure SaveFigure('Fig4') end %% (Part 5) Essential Functions that depend on the model function output = getR0 % Returns R0 given the parameter values % Find the formula for R0 on paper first, then code it here output =sqrt((mu1+omega)*beta1*beta2*gamma*delta/(d1*d2*d3*d4*mu1*mu3*(mu1+nu+omega))); end function output = getR0vaccination(mynu) % Returns R0 given the parameter values nu = mynu; CalculateAuxiliaryVariables; output =getR0; end function sols = getEquilibrium % Returns steady states of the dynamics % Find the formulas for the solutions on paper first, then code it % here. % The equilibrium will likely depend on R0, so we have to find that % first R0 = getR0; if R0<= 1 % input disease free ansS1=(mu1+omega)/(mu1*(mu1+nu+omega)); ansS2=(1/d2); ansV1=(nu/(mu1*(mu1+nu+omega))); ansI1=0; ansI2=0; ansM=0; ansP=0; else % input endemic equil. x = mu1*(omega+mu1+nu)/(omega+mu1); % x is 1/S_1^0 y = beta2*gamma*mu3; a1= x*d1*d2*d3*d4*alpha2*gamma.^2; a2= y*(delta*beta1*(d1-eta) + d1*d3*d4*x + d1*alpha1*delta*x); a3= x*d1*d2*d3*d4*mu3.^2 - delta*beta1*y; ansI1=(-a2+sqrt(a2.^2-4*a1*a3))/(2*a1); ansS1=(1-(d1-eta)*ansI1)/x; ansV1=nu/(omega+mu1)*ansS1; ansM = gamma/mu3*ansI1; ansS2= (1 + alpha2*ansM^2)./(beta2*ansM+d2*(1+alpha2*ansM^2)); ansI2 = beta2*ansM*ansS2/(d3*(1+alpha2*ansM^2)); ansP = delta*ansI2/d4; end % Now, with I done, solve for the rest of the compartments % Set the output (in the same order as the compartments are % ordered) sols = [ansS1, ansI1, ansV1, ansS2, ansI2, ansM, ansP]; end function output = getIncentiveFunction % Returns the incentive function, i.e. the cost of not-vaccination % minus the cost of vaccination % Get analytical solutions to see what the equilibria are sols = getEquilibrium; % Unpack to get the I as the second coordinate myP = sols(7); piNV = beta1*myP/(1+alpha1*myP)/(beta1*myP/(1+alpha1*myP)+mu1); output = piNV * CD - (omega/(omega+mu1) * piNV * CD +CV); end function output = getpNE % Returns Nash equilibrium, i.e. the root of the incentive function % Solve for pNE as a root of the incentive function output = getRoot(@getIncentiveFunction, 'nu', 0, 40, 0.00001); end function output = getINCne % Returns the incidence rate when everyone uses the nash % equilibrium % Get NE nuNE = getpNE; output = getINC(nuNE); end function output = getINC(mynu) % Returns incidence rate at vaccination rate mynu % Set vaccination rate to be NE nu = mynu; CalculateAuxiliaryVariables % Get equlibrium values sols = getEquilibrium; % P is the last coordinate myP = sols(end); myS = sols(1); %incidence rate is output= beta1*myP/(1+alpha1*myP) * myS/lambda1/(sum(sols(1:3))) * 10^5; end function output = getPRne % Returns the prevalence when everyone uses the nash % equilibrium % Get NE nuNE = getpNE; output = getPR(nuNE); end function output = getPR(mynu) % Returns prevalence at vaccination rate mynu % Set vaccination rate to be NE nu = mynu; CalculateAuxiliaryVariables % Get equlibrium values sols = getEquilibrium; %incidence rate is output= sols(2)/(sum(sols(1:3))); end function justTry %Set minimal error, firt to something huge minerr = 1000; for i=1:10000 beta1 = rand(1)*0.002+0.00; beta2 = rand(1)*10+10; alpha1 = rand(1)*0.05+0.0; M0 = rand(1)*10000; epsilon = rand(1)*10+1; err = getError([beta1, beta2, alpha1, M0, epsilon]); if err1 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((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';']) end % Calculate the auxiliary variables CalculateAuxiliaryVariables % Calculate Equilibrium Solutions ansols = getEquilibrium; % 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.I>0))) % if R0<1 but there is an endemic equilibrium (positive infected) disp('R0<1 but endemic equilibrium exists') displayParameters end if ((Reff>1) && (all(numsols.I<=0))) % if R0>1 but there is no endemic equilibrium (no positive infected) disp('R0>1 but endemic equilibrium does not exist') displayParameters end % Set the number of numerical solutions (initialize first) %n_of_numsols = []; %eval(['n_of_numsols = size(numsols.' Comps{1} ',1);']); % This finds the size of the solutions of the first compartment n_of_numsols = size(numsols.(Comps{1}),1); % This is a fix of above suggested by matlab % 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.I(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]; % % 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 between numerical 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 output = getnuHI(tolerance) % Returns the herd imunity value of p, i.e., where R0 is 1 % tolerance specifies how precise we want the pHI to be % it is set to 0.001 by default which means that pHI is good to 3 % decimal places. if nargin<1 % if no argument is given % Set the tolerance to 0.001 by default tolerance = 0.00001; end % The herd immunity happens where R0=1, i.e. the root of getR0-1=0 output = getRoot(@() getR0-1, 'nu', 0, 40, tolerance); % @() getR0-1 creates a handle to (an anonymous) function that returns R0-1 end function output = getCriticalValue(aa,bb, VarName) % Returns a critical value of Varname within interval [aa, bb] % i.e. value where R0 = 1 % Implicitly assumes that R0 is decreasing in the parameter VarName % but will work on increasing as well. May not work properly if teh % relationship is more complicated % Uses bisection method (recursively divides the intervals in half % until the interval is short enough and within our tolerance for % an error % aa should be smaller than bb (but this requirement is not % crucial especially if R0 is not decreasing) %THIS FUNCTION MAY BE OBSOLETE AND CAN LIKELY BE REPLACED BY GETROOT FUNCTION if abs(aa-bb)<10^(-3) %if we are already close enough together % return the midpoint output = (aa+bb)/2; else % we are still far away, so evaluate R0 % Store the current value of VarName %eval(['temp =' VarName ';']); % like temp = p; temp =eval(VarName); % Evaluate R0 at aa and bb and the midpoint eval([VarName '= aa ;']); % like p = aa; CalculateAuxiliaryVariables R0aa = getR0; eval([VarName '= bb ;']); % like p = bb; CalculateAuxiliaryVariables R0bb = getR0; eval([VarName '= (aa + bb)/2;']); % like p = (aa+bb)/2; CalculateAuxiliaryVariables R0mid = getR0; % Return VarName to its original value eval([VarName ' = temp ;']); % like p = temp; CalculateAuxiliaryVariables if (R0aa-1)*(R0mid-1)<0 % if R0 at aa has a different sign than R0 at midpoint output = getCriticalValue(aa,(aa+bb)/2, VarName); elseif (R0bb-1)*(R0mid-1)<0 % if R0 at bb has a different sign than R0 at midpoint output = getCriticalValue((aa+bb)/2,bb, VarName); else % should not happen, but if it does, it means that we cannot % use the method % In general, the function R0 is decreasing in our variable, so as long % as both are positive, it means that the critical values % is either the lower bound or the upper bound if R0aa<=1 % R0 already small enough at lower bound output = aa; elseif R0bb>=1 % R0 still too large at upper bound output = bb; else % This should really not happen but display warning % just in case disp('Something went wrong when calculating critical values') displayParameters end end end end function output = getRoot(myfunction, VarName, lb, ub, tol) % Returns a root of myfunction as a function of a variable Varname on interval [lb, ub] % the root is returned up to a tolerance value tol % % Implicitly assumes that myfunction is decreasing (but works in most % other cases too, especially if the first call is set up properly with myfunction % at aa having a different sign than at bb) % Uses bisection method (recursively divides the intervals in half % until the interval is short enough and within our tolerance for % an error % aa should be smaller than bb for decreasing functions but not % crucial if the sign of myfunction at aa is different than the % sign at bb % Declare the arguments and set their values if not provided if nargin <5 % if no tolerance is given % Set it as 10^-3 tol = 10^(-3); end if abs(lb-ub)=1 % myfunction still too large at upper bound output = ub; else % This should really not happen but display warning % just in case disp('Something went wrong when calculating root') %displayParameters % UNCOMMENT WHEN WE WANT TO SEE %WHAT PARAMETER VALUES CAUSED THE ISSUE SO THAT WE %CAN INVESTIGATE WHAT IS WRONG % Return the midpoint in the spirit of better to return something than nothing (VERY BAD PROGRAMING) output = (lb+ub)/2; end end end end function plotCriticalValues(aa,bb,VarName, DispName) % Plots how R0 depends on Varname within interval [aa, bb] % aa should be smaller than bb % DispName is the way VarName should be displayed on the plot % Store the current value of VarName temp =eval(VarName); % Set all values of VarName. Make sure we include the critical % value too critVarName = getCriticalValue(aa,bb, VarName); allValues = sort(unique([linspace(aa,bb,101) critVarName])); % Initialize myR0 myR0 = zeros(1, length(allValues)); % For every value in [aa,bb] for auxVar = 1:length(allValues) % Evaluate R0 at corresponding value of VarName eval([VarName '= allValues(auxVar) ;']); % like p = aa; CalculateAuxiliaryVariables myR0(auxVar) = getR0; end % Return VarName to its original value eval([VarName ' = temp ;']); % like p = temp; CalculateAuxiliaryVariables % Plot the figure figure % Dependence of R0 on VarName plot(allValues, myR0, 'k-') hold on % Horizontal line at R0=1 plot([allValues(1) allValues(end)], [1 1], 'k:') % Vertical line at the critical value plot([critVarName critVarName ], [0, 1], 'k:') % Add labels, ticks, and format and save xlabel(DispName) ylabel('Reproduction number') xticks(linspace(aa,bb,6)) FormatFigure SaveFigure(['CriticalValues_' VarName]) end function prcc=PRCClocal(values, Y, FileName, YName, Pars) % PRCC procedure % input: values: values of the parameters % and 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; if nargin == 3 % if the last parameter is not specified % Do sensititivity on all parameters Pars = Params; end % Get the number of parameters, discard the number of columns [numPars,~] = size(Pars); % % Get the parameter labels myLabels = Pars(:,2); % Preallocate PRCC values prcc = zeros(1,numPars); for i=1:numPars % 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 only the important coefficients %figure % Determine which ones are important importantPRCC = sortedPRCC(abs(sortedPRCC)>0.05); if isempty(importantPRCC) importantPRCC = sortedPRCC; end 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 = myLabels(I(abs(sortedPRCC)>0.05)); ax.TickLabelInterpreter = 'latex'; ylabel('Parameters', 'Interpreter', 'latex'); % Set the x axis between -1 and 1 xlim([-1 1]); xlabel(['PRCC for ' YName]); % Make the color gray h.FaceColor = [0.5 0.5 0.5]; %FormatFigure; %SaveFigure([FileName '_bar_Marino_important']) end function LHS_PRCC(RespFunction, RespName, FileName, Pars, Samples, tile_for_hist, tile_for_bars) % Does PRCC uncertainty and sensitivity analysis % 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 % Set the parameters if some are not provided if nargin < 5 % if Samples not provided % Specify how many times we want to generate the parameter % values Samples = 1000; end if nargin < 4 % if the Pars is not specified % Do sensititivity on all parameters Pars = Params; end disp('Doing sensitivity analysis based on Marino et al') % Get the number of parameters, discard the number of columns [numPars,~] = size(Pars); % % Get the parameter names myNames = Pars(:,1); % Initiate paramsUn = zeros(Samples, numPars); % This stores the values of the parameters allResp_Values = zeros(1,Samples); % This stores the values of the response function % Generate Latin Hypercube Sampling for aux1 = 1:Samples for aux = 1:numPars % Go through all the parameters % Evaluate the command that assigns the parameter with a given distribution with a parameter range % We need a column vector, but we will generate a row vector % first and transpose it during the assignment so that we can % use the eval command % Assign parameters at random but within bounds eval([myNames{aux} '=' num2str((Pars{aux,5}-Pars{aux,4})*rand+Pars{aux,4}) ';']) % Evaluate the command that assigns the parameter with a given name the proper base value %eval(['paramsUn(aux1, aux) = ' myNames{aux} ';']) paramsUn(aux1, aux) = eval(myNames{aux}); % Calculate auxiliary variables CalculateAuxiliaryVariables end % Evaluate the response variable allResp_Values(aux1) = RespFunction(); end % Plot the histogram to show uncertainty %figure nexttile(tile_for_hist) %allResp_Values histogram(allResp_Values, 'Normalization', 'probability'); xlabel(RespName) ylabel('Frequency') nexttile(tile_for_bars) %FormatFigure %SaveFigure([FileName '_hist_Marino']) % Get tha average value, do not consided NaN (not a number) values disp(['The average value of ' RespName ' is ' num2str(nanmean(allResp_Values))]); % Plot the bars with PRCC indices PRCClocal(paramsUn, allResp_Values', FileName, RespName, Pars) end function 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 %% (Part 7) Auxiliary Functions for easy plotting, etc function FormatSubFigure % Formats figures so that all figures are uniformly formatted h = gcf; % Get axes handle haxes=get(h, 'CurrentAxes'); % 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); end 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') saveas(h, ['Figures/' filename],'epsc') % 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