Corrected low_k_cutoff to be in line with paper definition

docs/src/development/implementation-details.md

@@ -42,10 +42,10 @@ If the solver tolerances need to be modified during adaptive tolerance solution
 
 ## [Removing Low-Rate Reactions](@id implementation_low_rate)
 
-When assembling a generated CRN into a system of ODEs to be integrated, some reactions can be purposefully left out if they have sufficiently low rates. This can be triggered with the `ODESimulationParams.low_k_cutoff` parameter by giving it a value other than `:none`. This calculates the maximum rate constant for every reaction, and reactions are then removed from the CRN if this maximum rate is below the cutoff. Setting the cutoff to `:auto` predicts a safe value, where maximum rates below this value will never change any species concentrations. This is checked by making sure that if a given reaction were to run at its maximum rate over the entire duration of the requested simulation (`ODESimulationParams.tspan[2]`, ``t_{\text{global}}``), any generated species concentration would still be less than the ODE solver's relative tolerance (`ODESimulationParams.reltol`):
+When assembling a generated CRN into a system of ODEs to be integrated, some reactions can be purposefully left out if they have sufficiently low rates. This can be triggered with the `ODESimulationParams.low_k_cutoff` parameter by giving it a value other than `:none`. This calculates the maximum rate constant for every reaction, and reactions are then removed from the CRN if this maximum rate is below the cutoff. Setting the cutoff to `:auto` predicts a safe value, where maximum rates below this value will never change any species concentrations. This is checked by making sure that if a given reaction were to run at its maximum rate over the entire duration of the requested simulation (`ODESimulationParams.tspan[2]`, ``t_{\text{global}}``) with maximum reactant concentrations set by `ODESimulationParams.low_k_maxconc` (``c^2_{\text{max}}``), any generated species concentration would still be less than the ODE solver's relative tolerance (`ODESimulationParams.reltol`):
 
 ```math
-\text{remove reaction } i \text{ if } t_{\text{global}} \cdot k_i^{\text{max}} < reltol
+\text{remove reaction } i \text{ if } t_{\text{global}} \cdot c^2_{\text{max}}k_i^{\text{max}} < reltol
 ```
 
 Maximum rates of reaction are calculated by enumeration over all permutations of maximum and minimum values of the conditions being calculated against. This should account for the vast majority of cases as rates are usually highest at extrema of experimental conditions.

docs/src/tutorials/ode-solution.md

@@ -235,11 +235,11 @@ This allows for controlling how often the ODE solver saves species concentration
 !!! warning "Saving with chunkwise time"
     When `solve_chunks=true`, this parameter behaves slightly differently. If left empty, species concentrations will be saved at the start and end of every chunk. If given a numeric value, this will be interpreted as a value within the local timescale `0.0:save_interval:solve_chunkstep`, so providing a value larger than `solve_chunkstep` will cause nothing to be saved. See the section on [Implemetation Details - Chunkwise Time](@ref implementation_chunkwise_time) for further details.
 
-#### `low_k_cutoff`
+#### `low_k_cutoff` and `low_k_maxconc`
 
-This parameter allows Kinetica to automatically remove reactions with low rate constants from a CRN, optimising CRN compilation and solution. It defaults to `low_k_cutoff=:auto`, which selects a conservative value for the minimum rate constant which ensures that only reactions which could never contribute to the current kinetic simulation are removed. 
+These parameters allows Kinetica to automatically remove reactions with low rate constants from a CRN, optimising CRN compilation and solution. They default to `low_k_cutoff=:auto` and `low_k_maxconc=2.0`, which selects conservative values that ensure that only reactions which could never contribute to the current kinetic simulation are removed. 
 
-If a numeric value is provided, this is used as the minimum rate constant below which reactions are removed from the CRN. If `low_k_cutoff=:none`, this behaviour is disabled and no reactions are removed from the CRN. See the section on [Implementation Details - Removing Low-Rate Reactions](@ref implementation_low_rate)
+If a numeric value is provided to `low_k_cutoff`, this is used as the minimum rate constant below which reactions are removed from the CRN. If `low_k_cutoff=:none`, this behaviour is disabled and no reactions are removed from the CRN. See the section on [Implementation Details - Removing Low-Rate Reactions](@ref implementation_low_rate)
 
 #### `allow_short_u0`
 

src/solving/params.jl

@@ -20,6 +20,7 @@ Contains fields for:
 * Whether to display progress bars - requires `TerminalLogger` initialisation (`progress=false`)
 * Time interval to interpolate solution data on (`save_interval=nothing`)
 * Cutoff below which reactions with low rate constants are removed from the network (`low_k_cutoff=:auto`)
+* Maximum species concentrationto multiply maximum reaction rates by un low rate cutoff (`low_k_maxconc=2.0`)
 * Whether to allow a vector `u0` to be shorter than the number of species in the network (`allow_short_u0=false`)
 """
 @kwdef mutable struct ODESimulationParams{tType, uType} <: AbstractSimulationParams
@@ -44,5 +45,6 @@ Contains fields for:
     # Optional network parameters
     save_interval::Union{tType, Nothing}=nothing
     low_k_cutoff::Union{uType, Symbol}=:auto
+    low_k_maxconc::uType = 2.0
     allow_short_u0::Bool=false
 end

src/solving/solve_utils.jl

@@ -202,6 +202,11 @@ If the cutoff is a numeric value, it is used directly. If it is
 reactions would not contribute to the network over the timespan
 of the simulation. If it is `:none`, does not apply a cutoff.
 
+Multiplies the calculated maximum rates of reaction by a theoretical
+maximum concentration that any reactants could attain through a 
+kinetic simulation, as set by `pars.low_k_maxconc`. This concentration
+multiplier is squared to emulate a bimolecular reaction.
+
 Returns the number of low-rate reactions that have been removed
 from the CRN.
 """
@@ -222,7 +227,7 @@ function apply_low_k_cutoff!(rd::RxData{iType, fType}, calc::cType,
     end
 
     # Calculate maximum rate constants.
-    max_rates = get_max_rates(conditions, calc)
+    max_rates = get_max_rates(conditions, calc) .* pars.low_k_maxconc^2
 
     # Find low rate reactions.
     low_rate_ids = Vector{iType}()