diff --git a/cara/montecarlo.py b/cara/montecarlo.py
index 4f8dc4e9..75c04771 100644
--- a/cara/montecarlo.py
+++ b/cara/montecarlo.py
@@ -818,7 +818,7 @@ def generate_boxplot(masked: bool = False, samples: int = 200000, qid: int = 100
 
 def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid: int = 100) -> None:
     """
-    Generates and displays a boxplot for comparing the qR-values in scenarios with various combinations of breathing
+    Generates and displays CDFs for comparing the qR-values in scenarios with various combinations of breathing
     rates and expiratory activities
     :param masked: A bool indicating whether or not the infected subject is wearing a mask
     :param samples: The number of samples to use for the monte carlo simulation
@@ -827,8 +827,8 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid:
     """
     fig, axs = plt.subplots(3, 1, sharex='all')
 
-    plt.suptitle("$F(qR|qID=$" + str(qid) + "$)$",fontsize=12, y=0.93)
-    fig.text(0.02, 0.5, '', va='center', rotation='vertical',fontsize=12)
+    plt.suptitle("$F(qR|qID=$" + str(qid) + "$)$",fontsize=14, y=0.93)
+    fig.text(0.02, 0.5, 'Cumulative Distribution Function', va='center', rotation='vertical',fontsize=14)
 
     scenarios = [MCInfectedPopulation(
         number=1,
@@ -844,7 +844,7 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid:
     right = max(np.max(qrs) for qrs in qr_values)
 
     # Colors can be changed here
-    colors = ['lightgrey', 'grey', 'black']
+    colors = ['skyblue', 'royalblue', 'navy']
 
     breathing_rates = ['Seated', 'Light activity', 'Heavy activity']
     activities = ['Breathing', 'Speaking', 'Shouting']
@@ -862,7 +862,7 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid:
         axs[i].grid(linestyle='--')
         axs[0].legend(handles=lines, loc='upper left')
 
-    plt.xlabel('qR ($q\;h^{-1}$)')
+    plt.xlabel('qR (q h$^{-1}$)', fontsize=12)
     tick_positions = np.arange(int(np.ceil(left)), int(np.ceil(right)), 2)
     plt.xticks(ticks=tick_positions, labels=['$\;10^{' + str(i) + '}$' for i in tick_positions])
 
@@ -1021,16 +1021,16 @@ def plot_concentration_curve(model: MCExposureModel):
         y[high_indexes] += lower_threshold + (upper_threshold - lower_threshold) * lower_scaling_factor + (x[high_indexes] - upper_threshold) * upper_scaling_factor
         return y
 
-    plt.xlabel('Time ($h$)', fontsize=12)
-    plt.ylabel('Concentration ($q\;m^{-3}$)', fontsize=12)
-    plt.title('Concentration of infectious quantum', fontsize=12)
+    plt.xlabel('Time (h)', fontsize=14)
+    plt.ylabel('Concentration (q m$^{-3}$)', fontsize=14)
+    plt.title('Concentration of infectious quantum', fontsize=14)
     plt.plot([start, stop], [upper_threshold, upper_threshold], linestyle='dotted', color='grey')
     plt.plot([start, stop], [lower_threshold, lower_threshold], linestyle='dotted', color='grey')
     plt.ylim(0, top)
 
     plt.yscale('function', functions=(forward, inverse))
 
-    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
+    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=10)
     plt.yticks([i for i in np.linspace(0, lower_threshold, 3)][:-1] +
                [i for i in np.linspace(lower_threshold, upper_threshold, 6)][:-1] +
                [i for i in np.linspace(upper_threshold, top, 3)])
@@ -1063,6 +1063,7 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label
 
     ax.legend(loc='upper left')
     ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2)
+    ax.spines["right"].set_visible(False)
 
     factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models]
     present_indexes = np.array([exp_models[0].exposed.person_present(t) for t in times])
@@ -1074,25 +1075,28 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label
            for c, factor in zip(modified_concentrations, factors)]
 
     plt.suptitle(title)
-    plt.xlabel("Time ($h$)", fontsize=12)
-    plt.ylabel("Quantum concentration ($q\;m^{-3}$)\nmean values of $C(t)$", fontsize=12)
+    plt.xlabel("Time ($h$)", fontsize=14)
+    plt.ylabel("Quantum concentration ($q\;m^{-3}$)\nmean values of $C(t)$", fontsize=14)
     if show_qd:
 
         ax1 = ax.twinx()
         for qd, label, color in zip(qds, labels, colors):
             ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted')
-        ax1.set_ylabel('Dose ($q$)', fontsize=12)
+        ax1.spines["right"].set_linestyle("--")
+        ax1.spines["right"].set_linestyle((0,(1,5)))
+        ax1.set_ylabel('Dose ($q$)', fontsize=14)
         ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2)
 
         ax2 = ax.twinx()
-        ax2.spines["right"].set_position(("axes", 1.2))
-        ax2.set_ylabel('Dose (RNA copies)\nmean values $qD$', fontsize=12)
+        ax2.spines["right"].set_position(("axes", 1.15))
+        ax2.spines["right"].set_linestyle((0,(1,5)))
+        ax2.set_ylabel('Dose (RNA copies)\nmean values $qD$', fontsize=14)
         ax2.set_ylim(tuple(y * exp_models[0].concentration_model.virus.qID for y in ax1.get_ylim()))
         #ax1.legend(loc='upper right')
 
     plt.tight_layout()
     plt.hlines([60], colors=['lightgrey'], linestyles=['dashed'], xmin=start, xmax=stop)
-    plt.text(7, 65, "$qID = 60$", color='lightgrey')
+    plt.text(7, 65, "$qID = 60$", color='grey')
     plt.show()