added partial mitigation for prep users to big model

Author: laura-mueller

notebooks/Nobs_Nreal_derivatives_bigmodel.py

@@ -21,7 +21,8 @@
 onlyplots = False # if True, only plots are generated, if False, first data is generated and then plots
 derivative = "dP" # dP: derivative with respect to PrEP adoption, dH: derivative with respect to risk awareness
 
-which_c = 8
+which_c = 9
+which_xi = 0
 
 sets_of_c = jnp.array([
     [31.0,  40.0,   60.0,  203.0],
@@ -32,7 +33,13 @@
     [2.0, 80.0, 53.3, 200.0],
     [2.0, 63.5, 100.0, 200.0], #6
     [75, 39, 18.5, 2], #7
-    [10.0, 50.0, 120.0, 190.0] #8
+    [10.0, 50.0, 120.0, 190.0], #8
+    [50.0, 50.0, 50.0, 50.0] #9
+])
+sets_of_xi = jnp.array([
+    [0.0, 0.0, 0.0, 0.0],
+    [0.8, 0.6, 0.4, 0.2], #1
+    [0.2, 0.4, 0.6, 0.8] #2
 ])
 
 
@@ -71,6 +78,7 @@ def calc_Nreal(H, P, lambdaP, betaSTI):
     args_mod["lambda_P"] = lambdaP
     args_mod["beta_STI"] = betaSTI
     args_mod["c"] = sets_of_c[which_c]
+    args_mod["xi"] = sets_of_xi[which_xi]  # partial mitigation for people on PrEP
 
     # run the model fro 80 years (long time to get to steady state)
     output = icomo.diffeqsolve(args = args_mod, ODE = bigmodel_STI.main_model, y0 = y0, ts_out = np.linspace(0, 365*80, 365*80+1), max_steps=365*80+1)
@@ -137,7 +145,7 @@ def calc_Nobs(H, P, lambdaP, betaSTI):
             # save stuff as npy files
 
             with open(
-                "../results/Nreal_Nobs_bigmodel_lambdap%g_betaSTI%g_dN%s_c%s.npy" %(lambdaP * 360, betaSTI, derivative, which_c),"wb",) as f:
+                "../results/Nreal_Nobs_bigmodel_lambdap%g_betaSTI%g_dN%s_c%s_xi%s.npy" %(lambdaP * 360, betaSTI, derivative, which_c, which_xi),"wb",) as f:
                 np.save(f, Ps)
                 np.save(f, Hs)
                 np.save(f, dNrealdP)
@@ -198,7 +206,7 @@ def discretize_cmaps(cmap, N):
     res = []
     for lambdaP in lambda_P_values:
         # load data
-        with open("../results/Nreal_Nobs_bigmodel_lambdap%g_betaSTI%g_dN%s_c%s.npy" %(lambdaP * 360, betaSTI, derivative, which_c),"rb") as f:
+        with open("../results/Nreal_Nobs_bigmodel_lambdap%g_betaSTI%g_dN%s_c%s_xi%s.npy" %(lambdaP * 360, betaSTI, derivative, which_c, which_xi),"rb") as f:
             Ps_plot = np.load(f)
             Hs_plot = np.load(f)
             dNrealdP_plot = np.load(f)
@@ -268,4 +276,4 @@ def discretize_cmaps(cmap, N):
 #plt.tight_layout()
 #plt.show()
 
-fig.savefig("../figures/final_figure_withinflux_bigmodel_DERIVATIVES_dN%s_c%s.pdf" %(derivative, which_c), format="pdf", bbox_inches="tight")
+fig.savefig("../figures/final_figure_withinflux_bigmodel_DERIVATIVES_dN%s_c%s_xi%s.pdf" %(derivative, which_c, which_xi), format="pdf", bbox_inches="tight")

src/tapm/HIVandSTI/bigmodel_STI.py

@@ -102,15 +102,16 @@ def duration2rate(x):
     [2.0, 80.0, 53.3, 200.0],
     [2.0, 63.5, 100.0, 200.0],
     [75, 39, 18.5, 2],
+    [50,50,50,50] #8
 ])
-c = sets_of_c[2] # 2, 5, 7
-sets_of_r = jnp.array([
-    [1.0,   1.0,    1.0,    1.0     ],
-    [1.5,   0.8,    0.3,    0.05    ],
-    [1.2,   1.1,    0.5,    0.113   ],
-    [0.5,   0.8,    2.0,    3.41    ]
+c = sets_of_c[8] # 2, 5, 7
+
+sets_of_xi = jnp.array([
+    [0.0, 0.0, 0.0, 0.0],
+    [0.8, 0.6, 0.4, 0.2], #1
+    [0.2, 0.4, 0.6, 0.8] #2
 ])
-r = sets_of_r[2] 
+xi = sets_of_xi[1]  # risk assimilation (partial mitigation for people on PrEP)
 
 #H = 5.0 # HIV hazard
 #P = 50.0 # PrEP fraction
@@ -121,7 +122,6 @@ def duration2rate(x):
     "mu": mu,
     "Omega": Omega,
     "c": c,
-    "r": r,
     "c_hiv": c_hiv,
     "h": h,
     "epsilon": epsilon,
@@ -148,7 +148,7 @@ def duration2rate(x):
     "gamma_STI": gamma_STI,
     "gammaT_STI": gammaT_STI,
     "beta_HIV": beta_HIV,
-
+    "xi": xi,
     #"H": H,
     #"P": P
     }
@@ -220,7 +220,8 @@ def foi_STI(y, args):
     """
 
     I_eff = y["Ia_STI"] + y["Is_STI"]
-    foi = args["beta_STI"] * (1 - m(args, y)*(1 - prep_fraction(y, args)) ) *contact_matrix(y, args) @ (I_eff/N_0)
+    # foi = args["beta_STI"] * (1 - m(args, y)*(1 - prep_fraction(y, args)) ) *contact_matrix(y, args) @ (I_eff/N_0) # without xi
+    foi = args["beta_STI"] * ((1 - m(args, y))*(1 - prep_fraction(y, args)) + (1-args["xi"]*m(args, y)) * prep_fraction(y, args)) *contact_matrix(y, args) @ (I_eff/N_0)
     return foi
 
 def contact_matrix(y, args):           
@@ -257,7 +258,7 @@ def prep_fraction(y, args):
     Returns:
         Value as float64.
     """
-    return args["P"]*args["r"]
+    return args["P"]
 
 
 def lambda_a(y, args):