diff --git a/cara/montecarlo.py b/cara/montecarlo.py
index 88553f39..947e9820 100644
--- a/cara/montecarlo.py
+++ b/cara/montecarlo.py
@@ -7,6 +7,29 @@ import typing
 import matplotlib.pyplot as plt
 from sklearn.neighbors import KernelDensity
 
+USE_SCOEH = False
+
+scoeh_vl_frequencies = ((1.880302953, 2.958422139, 3.308759599, 3.676921581, 4.036604757, 4.383770594,
+                         4.743136608, 5.094658141, 5.45613857, 5.812946142, 6.160090835, 6.518505362,
+                         6.866918705, 7.225333232, 7.574148314, 7.923640008, 8.283027166, 8.641758855,
+                         9.000448256, 9.35956054, 9.707720153, 10.06554264, 10.41435773, 10.76304594,
+                         11.12198907, 11.47118475),
+                        (0.001694915, 0.00720339, 0.027966102, 0.205932203, 0.213983051, 0.171186441,
+                         0.172881356, 0.217372881, 0.261440678, 0.211864407, 0.168644068, 0.151271186,
+                         0.133474576, 0.116101695, 0.106355932, 0.110169492, 0.112288136, 0.101271186,
+                         0.08940678, 0.086016949, 0.063135593, 0.033898305, 0.024152542, 0.011864407,
+                         0.005084746, 0.002966102))
+
+symptomatic_vl_frequencies = ((2.46032, 2.67431, 2.85434, 3.06155, 3.25856, 3.47256, 3.66957, 3.85979, 4.09927, 4.27081,
+                               4.47631, 4.66653, 4.87204, 5.10302, 5.27456, 5.46478, 5.6533, 5.88428, 6.07281, 6.30549,
+                               6.48552, 6.64856, 6.85407, 7.10373, 7.30075, 7.47229, 7.66081, 7.85782, 8.05653, 8.27053,
+                               8.48453, 8.65607, 8.90573, 9.06878, 9.27429, 9.473, 9.66152, 9.87552),
+                              (0.001206885, 0.007851618, 0.008078144, 0.01502491, 0.013258014, 0.018528495, 0.020053765,
+                               0.021896167, 0.022047184, 0.018604005, 0.01547796, 0.018075445, 0.021503523, 0.022349217,
+                               0.025097721, 0.032875078, 0.030594727, 0.032573045, 0.034717482, 0.034792991,
+                               0.033267721, 0.042887485, 0.036846816, 0.03876473, 0.045016819, 0.040063473, 0.04883754,
+                               0.043944602, 0.048142864, 0.041588741, 0.048762031, 0.027921732, 0.033871788,
+                               0.022122693, 0.016927718, 0.008833228, 0.00478598, 0.002807662))
 
 # The (k, lambda) parameters for the weibull-distributions corresponding to each quantity
 weibull_parameters = [((0.5951563631241763, 0.027071715346754264),          # emission_concentration
@@ -27,16 +50,7 @@ lognormal_parameters = ((0.10498338229297108, -0.6872121723362303),     # BR, se
                         (0.1894616357138137, 0.551771330362601),        # BR, moderate exercise
                         (0.21744554768657565, 1.1644665696723049))      # BR, heavy exercise
 
-log_viral_load_frequencies = ((1.880302953, 2.958422139, 3.308759599, 3.676921581, 4.036604757, 4.383770594,
-                               4.743136608, 5.094658141, 5.45613857, 5.812946142, 6.160090835, 6.518505362,
-                               6.866918705, 7.225333232, 7.574148314, 7.923640008, 8.283027166, 8.641758855,
-                               9.000448256, 9.35956054, 9.707720153, 10.06554264, 10.41435773, 10.76304594,
-                               11.12198907, 11.47118475),
-                              (0.001694915, 0.00720339, 0.027966102, 0.205932203, 0.213983051, 0.171186441,
-                               0.172881356, 0.217372881, 0.261440678, 0.211864407, 0.168644068, 0.151271186,
-                               0.133474576, 0.116101695, 0.106355932, 0.110169492, 0.112288136, 0.101271186,
-                               0.08940678, 0.086016949, 0.063135593, 0.033898305, 0.024152542, 0.011864407,
-                               0.005084746, 0.002966102))
+log_viral_load_frequencies = scoeh_vl_frequencies if USE_SCOEH else symptomatic_vl_frequencies
 
 
 def lognormal(csi: float, lamb: float, samples: int) -> np.ndarray: