init

Author: TorkelE

docs/src/assets/Project.toml

@@ -51,7 +51,7 @@ IncompleteLU = "0.2"
 JumpProcesses = "9.11"
 Latexify = "0.16"
 LinearSolve = "2.30"
-ModelingToolkit = "9.15"
+ModelingToolkit = "9.16.0"
 NonlinearSolve = "3.12"
 Optim = "1.9"
 Optimization = "3.25"
@@ -68,5 +68,5 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-StructuralIdentifiability = "0.5.7"
-Symbolics = "5.28"
+StructuralIdentifiability = "0.5.8"
+Symbolics = "5.30.1"

docs/src/faqs.md

@@ -59,7 +59,7 @@ plot(sol; idxs = [A, B])
 ```
 
 ## How to disable rescaling of reaction rates in rate laws?
-As explained in the [Reaction rate laws used in simulations](@ref) section, for
+As explained in the [Reaction rate laws used in simulations](@ref introduction_to_catalyst_ratelaws) section, for
 a reaction such as `k, 2X --> 0`, the generated rate law will rescale the rate
 constant, giving `k*X^2/2` instead of `k*X^2` for ODEs and `k*X*(X-1)/2` instead
 of `k*X*(X-1)` for jumps. This can be disabled when directly `convert`ing a
@@ -105,7 +105,7 @@ parametric_stoichiometry) section.
 
 ## How to set default values for initial conditions and parameters?
 How to set defaults when using the `@reaction_network` macro is described in
-more detail [here](@ref dsl_description_defaults). There are several ways to do
+more detail [here](@ref dsl_advanced_options_default_vals). There are several ways to do
 this. Using the DSL, one can use the `@species` and `@parameters` options:
 ```@example faq3
 using Catalyst

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -1,7 +1,7 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
 In this tutorial we provide an introduction to using Catalyst to specify
 chemical reaction networks, and then to solve ODE, jump, and SDE models
-generated from them. At the end we show what mathematical rate laws and
+generated from them[^1]. At the end we show what mathematical rate laws and
 transition rate functions (i.e. intensities or propensities) are generated by
 Catalyst for ODE, SDE and jump process models.
 
@@ -151,8 +151,7 @@ underlying problem.
     variable-based parameter mappings, `u₀symmap` and `psymmap`, while when
     directly passing `repressilator` we could use either those or the
     `Symbol`-based mappings, `u₀map` and `pmap`. `Symbol`-based mappings can
-    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref),
-    see the [Basic Syntax](@ref basic_examples) section for more details.
+    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref).
 
 
 !!! note

docs/src/inverse_problems/behaviour_optimisation.md

@@ -1,5 +1,5 @@
 # [Optimization for non-data fitting purposes](@id behaviour_optimisation)
-In previous tutorials we have described how to use [PEtab.jl](@ref petab_parameter_fitting) and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
+In previous tutorials we have described how to use PEtab.jl and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
 
 ## [Maximising the pulse amplitude of an incoherent feed forward loop](@id behaviour_optimisation_IFFL_example)
 Incoherent feedforward loops (network motifs where a single component both activates and deactivates a downstream component) are able to generate pulses in response to step inputs[^2]. In this tutorial we will consider such an incoherent feedforward loop, attempting to generate a system with as prominent a response pulse as possible.

docs/src/inverse_problems/global_sensitivity_analysis.md

@@ -1,6 +1,6 @@
 # [Global Sensitivity Analysis](@id global_sensitivity_analysis)
 *Global sensitivity analysis* (GSA) is used to study the sensitivity of a function's outputs with respect to its input[^1]. Within the context of chemical reaction network modelling it is primarily used for two purposes:
-- [When fitting a model's parameters to data](@ref petab_parameter_fitting), it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
+- When fitting a model's parameters to data, it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
 - [When measuring some system behaviour or property](@ref behaviour_optimisation), it can help determine which parameters influence that property. E.g. for a model of a biofuel-producing circuit in a synthetic organism, GSA could determine which system parameters have the largest impact on the total rate of biofuel production.
 
 GSA can be carried out using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package. This tutorial contains a brief introduction of how to use it for GSA on Catalyst models, with [GlobalSensitivity providing a more complete documentation](https://docs.sciml.ai/GlobalSensitivity/stable/).
@@ -52,7 +52,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0
 
 !!! note
     We should make a couple of notes about the example above:
-    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_logarithmic_scale), this is advantageous in the context of inverse problems such as this one.
+    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. This is advantageous in the context of inverse problems such as this one.
     - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set.
     - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`.
     - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -129,13 +129,13 @@ If we from previous knowledge know that $kD = 0.1$, and only want to fit the val
 fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2]))
 nothing # hide
 ```
-Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_remake)). Now we create our modified loss function:
+Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_problems_remake)). Now we create our modified loss function:
 ```@example diffeq_param_estim_1 
 loss_function_fixed_kD = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); prob_generator = fixed_p_prob_generator, maxiters=10000, verbose=false, save_idxs=4)
 nothing # hide
 ```
 
-We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$* and $kP$):
+We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$ and $kP$):
 ```@example diffeq_param_estim_1 
 optprob_fixed_kD = OptimizationProblem(loss_function_fixed_kD, [1.0, 1.0])
 optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead())
@@ -160,7 +160,7 @@ nothing # hide
 corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`).
 
 ## [Parameter fitting to multiple experiments](@id optimization_parameter_fitting_multiple_experiments)
-Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
+Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref ensemble_simulations). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
 
 ## [Optimisation solver options](@id optimization_parameter_fitting_solver_options)
 Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument:

docs/src/model_simulation/ode_simulation_performance.md

@@ -8,7 +8,7 @@ Generally, this short checklist provides a quick guide for dealing with ODE perf
 1. If performance is not critical, use [the default solver choice](@ref ode_simulation_performance_solvers) and do not worry further about the issue.
 2. If improved performance would be useful, read about solver selection (both in [this tutorial](@ref ode_simulation_performance_solvers) and [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)) and then try a few different solvers to find one with good performance.
 3. If you have a large ODE (approximately 100 variables or more), try the [various options for efficient Jacobian computation](@ref ode_simulation_performance_jacobian) (noting that some are non-trivial to use, and should only be investigated if truly required).
-4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref investigating) (with preference for the former, which is easier to use).
+4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref ode_simulation_performance_parallelisation) (with preference for the former, which is easier to use).
 
 ## [Regarding stiff and non-stiff problems and solvers](@id ode_simulation_performance_stiffness)
 Generally, ODE problems can be categorised into [*stiff ODEs* and *non-stiff ODEs*](https://en.wikipedia.org/wiki/Stiff_equation). This categorisation is important due to stiff ODEs requiring specialised solvers. A common cause of failure to simulate an ODE is the use of a non-stiff solver for a stiff problem. There is no exact way to determine whether a given ODE is stiff or not, however, systems with several different time scales (e.g. a CRN with both slow and fast reactions) typically generate stiff ODEs.

docs/src/model_simulation/simulation_introduction.md

@@ -275,7 +275,7 @@ nothing # hide
 ```
 If the `@default_noise_scaling` option is used, that term is only applied to reactions *without* `noise_scaling` metadata.
 
-While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the [`remake_reactionsystem`](@ref) function can be used to achieve a similar effect.
+While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the [`remake_reactionsystem`](@ref simulation_structure_interfacing_problems_remake) function can be used to achieve a similar effect.
 
 ## [Performing jump simulations using stochastic chemical kinetics](@id simulation_intro_jumps)
 

docs/src/model_simulation/simulation_plotting.md

@@ -58,7 +58,7 @@ can be used to plot `X` only. When only a single argument is given, the vector f
 plot(sol; idxs = brusselator.X + brusselator.Y)
 ```
 
-## [Multi-plot plots](@id simulation_plotting_options)
+## [Multi-plot plots](@id simulation_plotting_options_subplots)
 It is possible to save plots in variables. These can then be used as input to the `plot` command. Here, the plot command can be used to create plots containing multiple plots (by providing multiple inputs). E.g. here we plot the concentration of `X` and `Y` in separate subplots:
 ```@example simulation_plotting
 plt_X = plot(sol; idxs = [:X])
@@ -71,7 +71,7 @@ When working with subplots, the [`layout`](https://docs.juliaplots.org/latest/la
 plot(plt_X, plt_Y; layout = (2,1), size = (700,500))
 ```
 
-## [Saving plots](@id simulation_plotting_options)
+## [Saving plots](@id simulation_plotting_options_saving)
 Once a plot has been saved to a variable, the `savefig` function can be used to save it to a file. Here we save our Brusselator plot simulation (the first argument) to a file called "saved_plot.png" (the second argument):
 ```@example simulation_plotting
 plt = plot(sol)
@@ -80,7 +80,7 @@ rm("saved_plot.png") # hide
 ```
 The plot file type is automatically determined from the extension (if none is given, a .png file is created).
 
-## [Phase-space plots](@id simulation_plotting_options)
+## [Phase-space plots](@id simulation_plotting_options_phasespace)
 By default, simulations are plotted as species concentrations over time. However, [phase space](https://en.wikipedia.org/wiki/Phase_space#:~:text=In%20dynamical%20systems%20theory%20and,point%20in%20the%20phase%20space.) plots are also possible. This is done by designating the axis arguments using the `idxs` option, but providing them as a tuple. E.g. here we plot our simulation in $X-Y$ space:
 ```@example simulation_plotting
 plot(sol; idxs = (:X, :Y))

docs/src/model_simulation/simulation_structure_interfacing.md

@@ -162,7 +162,7 @@ get_S(oprob)
 ```
 
 ## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation)
-When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref rbasic_CRN_library_two_statesef) programmatically, and use its symbolic variables to check, and update, an initial condition:
+When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref basic_CRN_library_two_states) programmatically, and use its symbolic variables to check, and update, an initial condition:
 ```@example structure_indexing_symbolic_variables
 using Catalyst
 t = default_t()