Automating the derivation of sensitivity indices

I am currently reviewing some papers on mathematical epidemiology and trying to replicate certain results on those papers. I want to calculate the sensitivity indices of model parameters (of the paper https://pubmed.ncbi.nlm.nih.gov/32834653/) with respect to the basic reproduction number using the formula where p is a parameter and is the sensitivity index of p. I have already calculated the indices manually, but I would like to know how to automate this process using MATLAB. The formula for and the parameter values are given below. The sensitivity index values given in the paper are . Thank you in advance!

beta_s = 0.274;
alpha_a = 0.4775;
alpha_u = 0.695;
mu = 0.062;
q_i = 0.078;
gamma_a = 0.29;
eta_i = 0.009;
eta_u = 0.05;
R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

[Bas123]

Does anyone have any ideas on this question? Thanks!

Answer 1

[Torsten]

You could use

syms beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u
beta_s_num = 0.274;
alpha_a_num = 0.4775;
alpha_u_num = 0.695;
mu_num = 0.062;
q_i_num = 0.078;
gamma_a_num = 0.29;
eta_i_num = 0.009;
eta_u_num = 0.05;
R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu));
%S_beta_s as an example
dR0dbeta_s = diff(R0,beta_s)

dR0dbeta_s = 

S_beta_s = beta_s_num/subs(R0,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])*subs(dR0dbeta_s,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])

S_beta_s = 

1

[Bas123]

Hey, thanks a lot for your answer. It was really helpful. However, when I try to calculate the sensitivity index for the parameter alpha_a, it gives me a different value from the one given in the paper. I used the following code for that.

S_alpha_a = alpha_a_num/subs(R0,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])*subs(dR0dalpha_a,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])

[Torsten]

Maybe you use numerical values for the parameters different from those used in the paper.

syms beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u
beta_s_num = 0.274;
alpha_a_num = 0.4775;
alpha_u_num = 0.695;
mu_num = 0.062;
q_i_num = 0.078;
gamma_a_num = 0.29;
eta_i_num = 0.009;
eta_u_num = 0.05;
R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu));
%S_beta_s as an example
dR0dalpha_a = diff(R0,alpha_a)

dR0dalpha_a = 

S_alpha_a = alpha_a_num/subs(R0,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])*subs(dR0dalpha_a,[beta_s alpha_a alpha_u mu q_i gamma_a eta_i eta_u],[beta_s_num alpha_a_num alpha_u_num mu_num q_i_num gamma_a_num eta_i_num eta_u_num])

S_alpha_a = 

[Bas123]

Thank you very much for your answer! It helped me a lot. However, I can confirm that the values that I have given here are accurate. The thing is that when I used R and Mathematica to compute these indices, I was able to get the values given in the paper. So, I am wondering where the problem lies.

[Torsten]

And what is your MATHEMATICA code ?

[Bas123]

Hey, thank you! Here is the Mathematica code that I used.

f[betas_, alphaa_, alphau_, mu_, qi_, gammaa_, etau_] := (betas*alphaa)/(gammaa+mu) + (betas*alphau*gammaa*(1-qi))/((gammaa+mu)*(etau+mu))

R0 = f[0.274, 0.4775, 0.695, 0.062, 0.078, 0.29, 0.05]

Sbetas = (betas/R0)*D[f[betas, alphaa, alphau, mu, qi, gammaa, etau], betas]

SbetasValue = Sbetas/.{betas->0.274, alphaa->0.4775, alphau->0.695, mu->0.062, qi->0.078, gammaa->0.29, etau->0.05}

[Torsten]

I don't see a difference between MATLAB and MATHEMATICA code.

If the results are different, I'd compare intermediate results, e.g. R0 and dR0dbeta_s, in both programs.