diff --git a/use_cases/diversity-vs-nti-intensity.jl b/use_cases/diversity-vs-nti-intensity.jl
index e5ab494ff..195e67a7a 100644
--- a/use_cases/diversity-vs-nti-intensity.jl
+++ b/use_cases/diversity-vs-nti-intensity.jl
@@ -23,11 +23,9 @@ n_foodweb = 50 # Replicates of trophic backbones.
 # and not before.
 # To do so, we have lowered the `abstol` and `reltol` arguments of
 # `TerminateSteadyState`.
-# Moreover, we have set the `verbose` argument of the `ExtinctionCallback`
-# to remove printed information about species extinctions.
 tmax = 10_000
-callback = CallbackSet(ExtinctionCallback(1e-5, false), TerminateSteadyState(1e-8, 1e-6))
-verbose = false # Do not show '@info' message when species go extinct.
+verbose = false # Do not show '@info' messages during simulation run.
+callback = CallbackSet(ExtinctionCallback(1e-5, verbose), TerminateSteadyState(1e-8, 1e-6))
 
 # Set up non-trophic interactions.
 intensity_range_size = 5
@@ -42,59 +40,56 @@ interaction_names = keys(intensity_values_dict)
 n_interaction = length(interaction_names)
 
 # Main simulation loop.
-# Each thread writes in its own DataFrame. We merge them at the end of the loop.
-df_to_merge = DataFrame[]
-Threads.@threads for i in 1:n_foodweb # We parallelize computation if possible.
-    df_tmp = DataFrame(;
+# Each thread handle a single food web and writes in its own DataFrame.
+# Merge the DataFrames at the end of the loop.
+dfs = [DataFrame() for _ in 1:n_foodweb] # Fill the vector with empty DataFrames.
+Threads.@threads for i in 1:n_foodweb # Parallelize computation if possible.
+    df_thread = DataFrame(;
         foodweb_id = Int64[],
         interaction = Symbol[],
         intensity = Float64[],
         diversity = Int64[],
     )
-    # First we create a trophic backbone.
+    # First, create a trophic backbone.
     foodweb = FoodWeb(nichemodel, S; C, Z)
-    # Secondly we loop on the different non-trophic interactions we want to study.
+    # Secondly, loop on the different non-trophic interactions to study.
     for interaction in interaction_names
         intensity_values = intensity_values_dict[interaction]
-        # Third, we increase the non-trophic interaction intensity.
+        # Third, increase the non-trophic interaction intensity.
         for intensity in intensity_values
-            # We add the given non-trophic interaction to the trophic backbone
+            # Add the given non-trophic interaction to the trophic backbone
             # with the given intensity.
             net = MultiplexNetwork(foodweb; [interaction => (L = L_nti, I = intensity)]...)
             functional_response = ClassicResponse(net; c = i_intra)
             params = ModelParameters(net; functional_response)
-            solution = simulate(params, rand(S); tmax, callback, verbose)
-            # We save the final diversity at equilibrium.
-            push!(df_tmp, [i, interaction, intensity, richness(solution[end])])
+            solution = simulate(params, rand(S); tmax, callback)
+            # Save the final diversity at equilibrium.
+            push!(df_thread, [i, interaction, intensity, richness(solution[end])])
         end
         @info "Interaction $interaction (foodweb $i) done."
     end
-    push!(df_to_merge, df_tmp)
+    dfs[i] = df_thread
 end
-df = reduce(vcat, df_to_merge)
 @info "All simulations done."
+df = reduce(vcat, dfs)
 
-# First, for each trophic backbone we compute the relative diversity difference
-# compared to the reference that is given by the trophic network that does not
-# contain non-trophic interactions, i.e. intensity = 0.
+# Compute the relative diversity different for each network
+# compared to the reference network that does not contain non-trophic interactions,
+# i.e. non-trophic intensity is null.
 groups = groupby(df, [:foodweb_id, :interaction])
-relative_difference(x, x_ref) = (x - x_ref) / x_ref
 for sub_df in groups
-    reference_row = findfirst(x -> x == 0, sub_df.intensity) # Trophic interactions only.
+    reference_row = findfirst(==(0), sub_df.intensity) # Trophic interactions only.
     S_ref = sub_df[reference_row, :diversity]
-    transform!(
-        sub_df,
-        :diversity => ByRow(S -> relative_difference(S, S_ref)) => :delta_diversity,
-    )
+    relative_difference(S) = (S - S_ref) / S_ref
+    transform!(sub_df, :diversity => ByRow(relative_difference) => :delta_diversity)
 end
-df = combine(groups, :) # Secondly, we recombine the DataFrames.
-# Third, for each interaction and intensity combination we compute the mean
-# and the confidence interval of diversity variation.
-ic95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
+df = combine(groups, :) # Recombine DataFrames.
+# Compute the mean and the confidence interval of diversity variation for each interaction.
+ci95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
 df_processed = combine(
     groupby(df, [:interaction, :intensity]),
     :delta_diversity => mean,
-    :delta_diversity => ic95,
+    :delta_diversity => ci95,
 )
 
 # Plot diversity variation versus non-trophic interaction intensity.
@@ -115,9 +110,9 @@ for (idx, interaction) in zip(indices, interaction_names)
     )
     x = df_interaction.intensity
     y = df_interaction.delta_diversity_mean
-    ic = df_interaction.delta_diversity_ic95
+    ci = df_interaction.delta_diversity_ci95
     scatterlines!(x, y)
-    errorbars!(x, y, ic; whiskerwidth = 5)
+    errorbars!(x, y, ci; whiskerwidth = 5)
 end
 
 # Write x and y labels.
diff --git a/use_cases/paradox-enrichment.jl b/use_cases/paradox-enrichment.jl
new file mode 100644
index 000000000..3b2aef96f
--- /dev/null
+++ b/use_cases/paradox-enrichment.jl
@@ -0,0 +1,62 @@
+import AlgebraOfGraphics: set_aog_theme!
+
+using CairoMakie
+using DataFrames
+using EcologicalNetworksDynamics
+
+"""
+    biomass_extrema(solution, last)
+
+Compute biomass extrema for each species during the `last` time steps.
+"""
+function biomass_extrema(solution, last)
+    trajectories = extract_last_timesteps(solution; last, quiet = true)
+    S = size(trajectories, 1) # Row = species, column = time steps.
+    [(min = minimum(trajectories[i, :]), max = maximum(trajectories[i, :])) for i in 1:S]
+end
+
+foodweb = FoodWeb([0 0; 1 0]); # 2 eats 1
+functional_response = ClassicResponse(foodweb; aᵣ = 1, hₜ = 1, h = 1);
+K_values = LinRange(1, 10, 50)
+tmax = 1_000 # Simulation length.
+verbose = false # Do not show '@info' messages during the simulation.
+df = DataFrame(;
+    K = Float64[],
+    B_resource_min = Float64[],
+    B_resource_max = Float64[],
+    B_consumer_min = Float64[],
+    B_consumer_max = Float64[],
+)
+
+# Run simulations: for each carrying capacity we compute the equilibrium biomass.
+for K in K_values
+    environment = Environment(foodweb; K)
+    params = ModelParameters(foodweb; functional_response, environment)
+    B0 = rand(2) # Inital biomass.
+    solution = simulate(params, B0; tmax, verbose)
+    extrema = biomass_extrema(solution, "10%")
+    push!(df, [K, extrema[1].min, extrema[1].max, extrema[2].min, extrema[2].max])
+end
+
+# Plot the orbit diagram with Makie.
+set_aog_theme!() # AlgebraOfGraphics theme.
+c_r = :green # Resource color.
+c_c = :purple # Consumer color.
+c_v = :grey # Vertical lines color.
+fig = Figure()
+ax = Axis(fig[2, 1]; xlabel = "Carrying capacity, K", ylabel = "Equilibrium biomass")
+resource_line = scatterlines!(df.K, df.B_resource_min; color = c_r, markercolor = c_r)
+scatterlines!(df.K, df.B_resource_max; color = c_r, markercolor = c_r)
+consumer_line = scatterlines!(df.K, df.B_consumer_min; color = c_c, markercolor = c_c)
+scatterlines!(df.K, df.B_consumer_max; color = c_c, markercolor = c_c)
+K0 = 2.3
+v_line1 = vlines!(ax, [K0]; color = c_v)
+v_line2 = vlines!(ax, [1 + 2 * K0]; color = c_v, linestyle = :dashdot)
+Legend(
+    fig[1, 1],
+    [resource_line, consumer_line, v_line1, v_line2],
+    ["resource", "consumer", "K₀", "1+2K₀"];
+    orientation = :horizontal,
+    tellheight = true, # Adjust the height of the legend sub-figure.
+    tellwidth = false, # Do not adjust width of the orbit diagram.
+)
diff --git a/use_cases/persistence-vs-producer-competition.jl b/use_cases/persistence-vs-producer-competition.jl
index 308831ecb..5fab321ca 100644
--- a/use_cases/persistence-vs-producer-competition.jl
+++ b/use_cases/persistence-vs-producer-competition.jl
@@ -1,9 +1,10 @@
-using AlgebraOfGraphics
+import AlgebraOfGraphics: set_aog_theme!
+import ColorSchemes: get, viridis
+import Statistics: mean, std
+
 using CairoMakie
-using ColorSchemes
 using DataFrames
 using EcologicalNetworksDynamics
-using Statistics
 
 # Define global community parameters.
 S = 20 # Species richness.
@@ -16,11 +17,9 @@ n_replicates = 100 # Number of food web replicates for each parameter combinatio
 C_values = [0.05, 0.1, 0.2] # Connectance values.
 tol_C = 0.01 # Tolerance on connectance when generating foodweb with the nichemodel.
 tmax = 2_000 # Simulation length.
-verbose = false # No information printed during simulation.
+verbose = false # Do not show '@info' messages during simulation run.
 
 """
-    standardize_K(foodweb, K, αij)
-
 Standardize total carrying capacity `K` for the number of producers in the `foodweb`
 and interspecific competition among producers `αij`.
 """
@@ -30,10 +29,10 @@ function standardize_K(foodweb, K, αij)
 end
 
 # Main simulation loop.
-# Each thread writes in its own DataFrame. We merge them at the end of the loop.
-df_to_merge = DataFrame[]
+# Each thread writes in its own DataFrame. Merge them at the end of the loop.
+dfs = [DataFrame() for _ in 1:length(C_values)] # Fill the vector with empty DataFrames.
 Threads.@threads for C in C_values # Parallelize on connctance values.
-    df_tmp = DataFrame(; C = Float64[], αij = Float64[], persistence = Float64[])
+    df_thread = DataFrame(; C = Float64[], αij = Float64[], persistence = Float64[])
     for i in 1:n_replicates
         foodweb = FoodWeb(nichemodel, S; Z, C, tol_C)
         for αij in αij_values
@@ -45,22 +44,22 @@ Threads.@threads for C in C_values # Parallelize on connctance values.
             # Measure species persistence i.e. the number of species that have
             # a biomass above 0 (`threshold`) at the last timestep (`last`).
             persistence = species_persistence(solution; threshold = 0, last = 1)
-            push!(df_tmp, [C, αij, persistence])
+            push!(df_thread, [C, αij, persistence])
             @info "C = $C, foodweb = $i, αij = $αij: done."
         end
     end
-    push!(df_to_merge, df_tmp)
+    push!(dfs, df_thread)
 end
-df = reduce(vcat, df_to_merge)
 @info "All simulations done."
+df = reduce(vcat, dfs)
 
-# We compute the mean persistence and the 95% confidence interval (`ic95`)
+# Compute the mean persistence and the 95% confidence interval (`ci95`)
 # for each (C, αij) combination.
 groups = groupby(df, [:C, :αij])
 df_processed = combine(
     groups,
     :persistence => mean,
-    :persistence => (x -> 1.96 * std(x) / sqrt(length(x))) => :ic95,
+    :persistence => (x -> 1.96 * std(x) / sqrt(length(x))) => :ci95,
 )
 
 # Plot mean species persistence with its confidence interval
@@ -73,17 +72,17 @@ ax = Axis(
     ylabel = "Species persistence",
 )
 curves = []
-colors = [get(ColorSchemes.viridis, val) for val in LinRange(0, 1, length(C_values))]
+colors = [get(viridis, val) for val in LinRange(0, 1, length(C_values))]
 for (C, color) in zip(C_values, colors)
     df_extract = df_processed[df_processed.C.==C, :]
     x = df_extract.αij
     y = df_extract.persistence_mean
-    ic = df_extract.ic95
+    ci = df_extract.ci95
     sl = scatterlines!(x, y; color = color, markercolor = color)
-    errorbars!(x, y, ic; color, whiskerwidth = 5)
+    errorbars!(x, y, ci; color, whiskerwidth = 5)
     push!(curves, sl)
 end
-Legend( # Put legend on top of the figure.
+Legend(
     fig[1, 1],
     curves,
     ["C = $C" for C in C_values];
diff --git a/use_cases/tritrophic-biomass-vs-temperature.jl b/use_cases/tritrophic-biomass-vs-temperature.jl
index d29091f95..0ea133071 100644
--- a/use_cases/tritrophic-biomass-vs-temperature.jl
+++ b/use_cases/tritrophic-biomass-vs-temperature.jl
@@ -2,8 +2,8 @@ import AlgebraOfGraphics: set_aog_theme!
 import ColorSchemes: tol_light
 
 using CairoMakie
-using EcologicalNetworksDynamics
 using DataFrames
+using EcologicalNetworksDynamics
 
 # Create tri-tropic network.
 tritrophic = FoodWeb([0 0 0; 1 0 0; 0 1 0])
@@ -20,6 +20,7 @@ T_values = values(T0:1:T40)
 
 # DataFrame to store the simulations outputs.
 df = DataFrame(; T = Float64[], trophic_level = Int64[], Beq = Float64[])
+verbose = false # Do not show '@info' messages during simulation run.
 
 # Run simulations for each temperature across gradient.
 for T in T_values
@@ -27,10 +28,10 @@ for T in T_values
     # Simulate biomass dynamics with modified parameters.
     B0 = params.environment.K[1] / 8 # Inital biomass.
     # `set_temperature!()` expresses rates in seconds,
-    # thus we also simulation time in seconds.
+    # then simulation time is expressed in seconds as well.
     # Simulation length and time steps are taken from Binzer et al., 2012.
     tmax = 315_360_000_000 # From Binzer et al., 2012.
-    callback = ExtinctionCallback(1e-6, false) # Remove TerminateSteadyState callback.
+    callback = ExtinctionCallback(1e-6, verbose) # Remove TerminateSteadyState callback.
     solution = simulate(params, [B0]; tmax, callback)
     Beq_vec = solution[end]
     for (Beq, tlvl) in zip(Beq_vec, trophic_lvl)
@@ -51,7 +52,7 @@ for (sub_df, color, markersize) in zip(groupby(df, :trophic_level), colors, mark
     sl = scatterlines!(T_celsius, sub_df.Beq; color, markercolor, markersize)
     push!(curves, sl)
 end
-Legend( # Put legend on top of the figure.
+Legend(
     fig[1, 1],
     curves,
     ["producer", "intermediate\nconsumer", "top predator"];
diff --git a/use_cases/trophic-classes-vs-mortality.jl b/use_cases/trophic-classes-vs-mortality.jl
index 8f112cc5e..ed4822c3c 100644
--- a/use_cases/trophic-classes-vs-mortality.jl
+++ b/use_cases/trophic-classes-vs-mortality.jl
@@ -1,8 +1,8 @@
-import ColorSchemes: tol_light # Color gradient for plot.
+import AlgebraOfGraphics: set_aog_theme!
+import ColorSchemes: tol_light
 import DiffEqCallbacks: CallbackSet, TerminateSteadyState
 import Statistics: mean, std
 
-using AlgebraOfGraphics
 using CairoMakie
 using DataFrames
 using DifferentialEquations
@@ -26,16 +26,15 @@ check_cycle = true # Check that generated food webs do not contain cycle(s).
 # and not before.
 # To do so, we have lowered the `abstol` and `reltol` arguments of
 # `TerminateSteadyState`.
-# Moreover, we have set the `verbose` argument of the `ExtinctionCallback`
-# to remove printed information about species extinctions.
-callback = CallbackSet(ExtinctionCallback(1e-5, false), TerminateSteadyState(1e-8, 1e-6))
+verbose = false # Do not show '@info' messages during simulation run.
+callback = CallbackSet(ExtinctionCallback(1e-5, verbose), TerminateSteadyState(1e-8, 1e-6))
 tmax = 10_000
 
 # Main simulation loop.
-# Each thread writes in its own DataFrame. We merge them at the end of the loop.
-df_to_merge = DataFrame[]
+# Each thread writes in its own DataFrame. Merge them at the end of the loop.
+dfs = [DataFrame() for _ in 1:n_foodweb] # Fill the vector with empty DataFrames.
 Threads.@threads for i in 1:n_foodweb # Parallelize computation when possible.
-    df_tmp = DataFrame(; # To store simulation data.
+    df_thread = DataFrame(; # To store simulation data.
         foodweb_id = Int64[],
         d0 = Float64[],
         n_producers = Int64[],
@@ -46,10 +45,11 @@ Threads.@threads for i in 1:n_foodweb # Parallelize computation when possible.
     # that are kept constant. We check that there is no cycle so that the trophic levels
     # and thus the trophic classes are well defined.
     foodweb = FoodWeb(nichemodel, S; C, Z, check_cycle)
-    functional_response = ClassicResponse(foodweb; c = i_intra) # Set intraspecific predator interference.
+    # Set intraspecific predator interference.
+    functional_response = ClassicResponse(foodweb; c = i_intra)
     B0 = rand(S) # Initial biomass associated with the food web.
     classes = nothing
-    # Then we loop on the mortality intensity that we increase at each iteration.
+    # Loop on the mortality intensity that is incremented at each iteration.
     for (j, d0) in enumerate(d0_values)
         # Scale default mortality rates by `d0`.
         d = d0 .* allometric_rate(foodweb, DefaultMortalityParams())
@@ -63,56 +63,52 @@ Threads.@threads for i in 1:n_foodweb # Parallelize computation when possible.
         surviving_sp = setdiff(1:richness(foodweb), extinct_sp)
         # Species classes are given by the the dynamics with zero mortality rate.
         d0 == 0 && (classes = trophic_classes(remove_species(foodweb, extinct_sp)))
-        # Then record the surviving species in each class
+        # Record the surviving species in each class
         # after running the dynamic with an increased mortality.
         n_per_class = [length(intersect(class, surviving_sp)) for class in classes]
-        push!(df_tmp, vcat(i, d0, n_per_class))
+        push!(df_thread, vcat(i, d0, n_per_class))
         @info "foodweb $i, d0 = $d0: done."
     end
-    push!(df_to_merge, df_tmp)
+    push!(dfs, df_thread)
 end
-df = reduce(vcat, df_to_merge)
 @info "All simulations done."
+df = reduce(vcat, dfs)
 
 # Process data.
-# First, we group the data per food web to compute for each food web
-# the relative difference in each trophic class compared to the reference (d0 = 0).
+# Group data per food web to compute the relative difference in diversity
+# for each trophic class of each food web,
+# compared to the reference given by d0 = 0.
 relative_difference(x, x_ref) = x_ref == 0 ? 0 : (x_ref - x) / x_ref
 groups = groupby(df, :foodweb_id)
 for sub_df in groups
-    reference_row = findfirst(x -> x == 0, sub_df.d0)
+    reference_row = findfirst(==(0), sub_df.d0)
+    ref_prod = sub_df[reference_row, :n_producers]
+    ref_cons = sub_df[reference_row, :n_intermediate_consumers]
+    ref_pred = sub_df[reference_row, :n_top_predators]
     transform!(
         sub_df,
         :d0,
-        :n_producers =>
-            ByRow(x -> relative_difference(x, sub_df[reference_row, :n_producers])) =>
-                :delta_producers,
+        :n_producers => ByRow(x -> relative_difference(x, ref_prod)) => :delta_producers,
         :n_intermediate_consumers =>
-            ByRow(
-                x -> relative_difference(
-                    x,
-                    sub_df[reference_row, :n_intermediate_consumers],
-                ),
-            ) => :delta_intermediate_consumers,
+            ByRow(x -> relative_difference(x, ref_cons)) => :delta_intermediate_consumers,
         :n_top_predators =>
-            ByRow(x -> relative_difference(x, sub_df[reference_row, :n_top_predators])) =>
-                :delta_top_predators,
+            ByRow(x -> relative_difference(x, ref_pred)) => :delta_top_predators,
     )
 end
-# Secondly, we recombine the DataFrames with the new column corresponding to the class
-# differences.
+
+# Recombine the DataFrames with the new column corresponding to the class differences.
 df = combine(groups, :)
-# Third, we group them per mortality rate (d0), and for each mortality we compute
-# the mean and confidence interval in each class difference.
-ic95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
+# Group DataFrames per mortality rate (d0), and compute the mean and the confidence
+# interval in each trophic class for each mortality value.
+ci95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
 df_processed = combine(
     groupby(df, :d0),
     :delta_producers => mean,
-    :delta_producers => ic95,
+    :delta_producers => ci95,
     :delta_intermediate_consumers => mean,
-    :delta_intermediate_consumers => ic95,
+    :delta_intermediate_consumers => ci95,
     :delta_top_predators => mean,
-    :delta_top_predators => ic95,
+    :delta_top_predators => ci95,
 )
 
 # Plot relative loss in each trophic class versus the mortality rate intensity.
@@ -124,14 +120,14 @@ x = df_processed[:, :d0]
 colors = tol_light[1:n_class]
 for (class, color) in zip(trophic_classes(), colors)
     col_mean = join(["delta", class, "mean"], "_")
-    col_ic = join(["delta", class, "ic95"], "_")
+    col_ci = join(["delta", class, "ci95"], "_")
     y = df_processed[:, col_mean]
-    ic = df_processed[:, col_ic]
+    ci = df_processed[:, col_ci]
     sl = scatterlines!(x, y; color, markercolor = color)
-    errorbars!(x, y, ic; color, whiskerwidth = 5)
+    errorbars!(x, y, ci; color, whiskerwidth = 5)
     push!(curves, sl)
 end
-Legend( # Put legend on top of the figure.
+Legend(
     fig[1, 1],
     curves,
     ["producers", "intermediate\nconsumers", "top predators"];