diff --git a/Short-Range_paper/results_paper.py b/Short-Range_paper/results_paper.py
index 4f6a6c2a8815d6bc5545883874ad3c1e187e8a82..685a69fc183774e11257d230aae70e49230be529 100644
--- a/Short-Range_paper/results_paper.py
+++ b/Short-Range_paper/results_paper.py
@@ -43,6 +43,8 @@ from model_scenarios_paper import *
 
 ############# FIGURES ############
 
+# Plots with viral loads and emission rates + statistical data ############
+# present_vl_er_histograms(physical_activity='Standing', expiratory_activities=['Breathing', 'Speaking', 'Shouting'])
 
 # Compare covid_overal_vl_data vs. symptomatic_vl_frequencies distributions
 # compare_vl_hist()
diff --git a/Short-Range_paper/scripts_paper.py b/Short-Range_paper/scripts_paper.py
index 628bf6e48bd9cc677f323a2ff2dc4f49035d46a0..e059e5d5b358167556f4f2d75805afc0a58460e3 100644
--- a/Short-Range_paper/scripts_paper.py
+++ b/Short-Range_paper/scripts_paper.py
@@ -18,6 +18,7 @@ import matplotlib.pyplot as plt
 from matplotlib.lines import Line2D
 import matplotlib as mpl
 from matplotlib.colors import rgb2hex
+
 import typing
 import dataclasses
 
@@ -865,8 +866,15 @@ def statistics_for_random_variable(physical_activity: bool = False,
     return None
 
 
-def short_range_emission_rate_values(physical_activity: str, expiratory_activity: str) -> None:
-    SAMPLE_SIZE = 250_000 if expiratory_activity == "Breathing" else 1_000_000
+def short_range_emission_rate_values(physical_activity: str, 
+                                     expiratory_activity: str,
+                                     virus_dist: typing.Optional[models.Virus] = None,
+                                     ):
+    
+    if virus_dist:
+        SAMPLE_SIZE = len(virus_dist.viral_load_in_sputum)
+    else:
+        SAMPLE_SIZE = 250_000 if expiratory_activity == "Breathing" else 1_000_000
 
     # Activity distribution
     activity_dist: models.Activity = activity_distributions(data_registry)[physical_activity].build_model(SAMPLE_SIZE)
@@ -878,7 +886,8 @@ def short_range_emission_rate_values(physical_activity: str, expiratory_activity
     short_range_distances_dist: models._VectorisedFloat = short_range_distances(data_registry).generate_samples(SAMPLE_SIZE)
 
     # Virus distribution
-    virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE)
+    if not virus_dist:
+        virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE)
 
     # Infected population
     infected: models.InfectedPopulation = models.InfectedPopulation(
@@ -908,4 +917,77 @@ def short_range_emission_rate_values(physical_activity: str, expiratory_activity
     print('\n<<<<<<<<<<< ' + expiratory_activity + ' model statistics >>>>>>>>>>>')
     print_er_info(emission_rate, log10_emission_rate)
 
-    return None
+    return emission_rate, log10_emission_rate
+
+def present_vl_er_histograms(physical_activity: str, expiratory_activities: typing.List) -> None:
+    # Virus distribution
+    virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE)
+
+    emission_rate_dict: typing.Dict = {}
+    log10_emission_rate_dict: typing.Dict = {}
+    for exp_act in expiratory_activities:
+        emission_rate, log10_emission_rate = short_range_emission_rate_values(physical_activity, exp_act, virus_dist)
+        emission_rate_dict[exp_act] = emission_rate
+        log10_emission_rate_dict[exp_act] = log10_emission_rate
+
+    # Viral load
+    viral_loads = models._VectorisedFloat = [np.log10(vl) for vl in virus_dist.viral_load_in_sputum]
+
+    # Plot material
+    fig, axs = plt.subplots(1, 2, sharex=False, sharey=False)
+    fig.set_figheight(5)
+    fig.set_figwidth(10)
+    plt.tight_layout()
+    plt.subplots_adjust(wspace=0.2)
+    plt.subplots_adjust(top=0.9)
+    plt.subplots_adjust(bottom=0.15)
+
+    axs[0].hist(viral_loads, bins = 300, color='lightgrey')
+    axs[0].set_xlabel('vl$_{\mathrm{in}}$ (log$_{10}$ RNA copies mL$^{-1}$)')
+    
+    # Viral load mean value
+    mean = np.mean(viral_loads)
+    axs[0].vlines(x=(mean), ymin=0, ymax=axs[0].get_ylim()[1], colors=('grey'), linestyles=('dashed'))
+
+    # Plot colors
+    hist_colors = ['lightsteelblue', 'wheat', 'darkseagreen']
+    hist_vlines_colors = ['cornflowerblue', 'goldenrod', 'olivedrab']
+    
+    # Plot the histograms
+    log10_emission_rate_dict_mean = {}
+    for exp_act_key, hist_color in zip(log10_emission_rate_dict, hist_colors):
+        exp_act_value = log10_emission_rate_dict[exp_act_key]
+        axs[1].hist(exp_act_value, bins=300, color=hist_color)
+        # Save the mean value
+        log10_emission_rate_dict_mean[exp_act_key] = np.mean(exp_act_value)
+
+    # Max value for y axis
+    ymax: float = axs[1].get_ylim()[1]
+
+    # Plot the vlines
+    for exp_act_key, vline_color in zip(log10_emission_rate_dict_mean, hist_vlines_colors):
+        exp_act_mean = log10_emission_rate_dict_mean[exp_act_key]
+        axs[1].vlines(x=exp_act_mean, ymin=0, ymax=ymax, colors=vline_color, alpha=0.75, linestyles='dashed')
+
+    axs[1].set_xlabel('vR (log$_{10}$ virions h$^{-1}$)')
+
+    from matplotlib.patches import Patch
+    import matplotlib.lines as mlines
+
+    hist_colors_legend: typing.List = ['lightgrey'] + hist_colors
+    hist_str_legend: typing.List = ['Viral Load'] + expiratory_activities
+
+    labels = [Patch(facecolor=[], edgecolor=[], color=color, label=label)
+             for color, label in zip(hist_colors_legend, hist_str_legend)]
+    labels.append(mlines.Line2D([], [], color='black',
+                            marker='', linestyle='dashed', label='Mean'))
+
+    for x in (0, 1):
+        axs[x].set_yticklabels([])
+        axs[x].set_yticks([])
+
+    plt.legend(handles=labels, loc='upper left', bbox_to_anchor=(1, 1))
+    plt.tight_layout()
+    plt.show()
+
+    return None
\ No newline at end of file