cleanup

Author: vyudu

docs/src/index.md

@@ -92,7 +92,7 @@ Pkg.add("Catalyst")
 
 Many Catalyst features require the installation of additional packages. E.g. for ODE-solving and simulation plotting
 ```julia
-Pkg.add("OrdinaryDiffEqTsit5")
+Pkg.add("OrdinaryDiffEqDefault")
 Pkg.add("Plots")
 ```
 is also needed.
@@ -124,7 +124,7 @@ an ordinary differential equation.
 
 ```@example home_simple_example
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEqTsit5, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 
 # Create model.
 model = @reaction_network begin
@@ -140,7 +140,7 @@ ps = [:kB => 0.01, :kD => 0.1, :kP => 0.1]
 ode = ODEProblem(model, u0, tspan, ps)
 
 # Simulate ODE and plot results.
-sol = solve(ode, Tsit5())
+sol = solve(ode)
 plot(sol; lw = 5)
 ```
 

docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md

@@ -55,15 +55,15 @@ To import a Julia package into a session, you can use the `using PackageName` co
 using Pkg
 Pkg.add("Catalyst")
 ```
-Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we also wish to install the needed sub-libraries of `OrdinaryDiffEqTsit5` and `Plots` packages (for numeric simulation of models, and plotting, respectively). We will import the default recommende dsolver, `Tsit5()`, which is in the `OrdinaryDiffEqTsit5` sub-library. A full list of `OrdinaryDiffEq` solver sublibraries can be found on the sidebar of [this page](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
+Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we also wish to install the needed sub-libraries of `OrdinaryDiffEq` and `Plots` packages (for numeric simulation of models, and plotting, respectively). We will import the default recommended solver from the `OrdinaryDiffEqDefault` sub-library. A full list of `OrdinaryDiffEq` solver sublibraries can be found on the sidebar of [this page](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 ```julia
-Pkg.add("OrdinaryDiffEqTsit5")
+Pkg.add("OrdinaryDiffEqDefault")
 Pkg.add("Plots")
 ```
 Once a package has been installed through the `Pkg.add` command, this command does not have to be repeated if we restart our Julia session. We can now import all three packages into our current session with:
 ```@example ex2
 using Catalyst
-using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqDefault
 using Plots
 ```
 Here, if we restart Julia, these `using` commands *must be rerun*.

docs/src/inverse_problems/global_sensitivity_analysis.md

@@ -22,7 +22,7 @@ end
 ```
 We will study the peak number of infected cases's ($max(I(t))$) sensitivity to the system's three parameters. We create a function which simulates the system from a given initial condition and measures this property:
 ```@example gsa_1
-using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqDefault
 
 u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0]
 p_dummy = [:β => 0.0, :a => 0.0, :γ => 0.0]
@@ -31,7 +31,7 @@ oprob_base = ODEProblem(seir_model, u0, (0.0, 10000.0), p_dummy)
 function peak_cases(p)
     ps = [:β => p[1], :a => p[2], :γ => p[3]]
     oprob = remake(oprob_base; p = ps)
-    sol = solve(oprob, Tsit5(); maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = false)
     SciMLBase.successful_retcode(sol) || return Inf
     return maximum(sol[:I])
 end

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -21,10 +21,10 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
 ps_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
 
 # Generate synthetic data.
-using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqDefault
 oprob_true = ODEProblem(rn, u0, (0.0, 10.0), ps_true)
-true_sol = solve(oprob_true, Tsit5())
-data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
+true_sol = solve(oprob_true)
+data_sol = solve(oprob_true; saveat=1.0)
 data_ts = data_sol.t[2:end]
 data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
 

docs/src/model_creation/conservation_laws.md

@@ -10,11 +10,11 @@ end
 ```
 If we simulate it, we note that while the concentrations of $X₁$ and $X₂$ change throughout the simulation, the total concentration of $X$ ($= X₁ + X₂$) is constant:
 ```@example conservation_laws
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol; idxs = [rs.X₁, rs.X₂, rs.X₁ + rs.X₂], label = ["X₁" "X₂" "X₁ + X₂ (a conserved quantity)"])
 ```
 This makes sense, as while $X$ is converted between two different forms ($X₁$ and $X₂$), it is neither produced nor degraded. That is, it is a *conserved quantity*. Next, if we consider the ODE that our model generates:
@@ -47,15 +47,15 @@ Here, Catalyst encodes all conserved quantities in a single, [vector-valued](@re
 
 Practically, the `remove_conserved = true` argument can be provided when a `ReactionSystem` is converted to an `ODEProblem`:
 ```@example conservation_laws
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps; remove_conserved = true)
 nothing # hide
 ```
 Here, while `Γ[1]` becomes a parameter of the new system, it has a [default value](@ref dsl_advanced_options_default_vals) equal to the corresponding conservation law. Hence, its value is computed from the initial condition `[:X₁ => 80.0, :X₂ => 20.0]`, and does not need to be provided in the parameter vector. Next, we can simulate and plot our model using normal syntax:
 ```@example conservation_laws
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 !!! note

docs/src/model_creation/dsl_advanced.md

@@ -100,20 +100,20 @@ end
 ```
 Next, if we simulate the model, we do not need to provide values for species or parameters that have default values. In this case all have default values, so both `u0` and `ps` can be empty vectors:
 ```@example dsl_advanced_defaults
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = []
 tspan = (0.0, 10.0)
 p = []
 oprob = ODEProblem(rn, u0, tspan, p)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 It is still possible to provide values for some (or all) initial conditions/parameters in `u0`/`ps` (in which case these overrides the default values):
 ```@example dsl_advanced_defaults
 u0 = [:X => 4.0]
 p = [:d => 0.5]
 oprob = ODEProblem(rn, u0, tspan, p)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 It is also possible to declare a model with default values for only some initial conditions/parameters:
@@ -127,7 +127,7 @@ end
 tspan = (0.0, 10.0)
 p = [:p => 1.0, :d => 0.2]
 oprob = ODEProblem(rn, u0, tspan, p)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 
@@ -142,6 +142,7 @@ end
 ```
 Please note that as the parameter `X₀` does not occur as part of any reactions, Catalyst's DSL cannot infer whether it is a species or a parameter. This must hence be explicitly declared. We can now simulate our model while providing `X`'s value through the `X₀` parameter:
 ```@example dsl_advanced_defaults
+using OrdinaryDiffEqTsit5
 u0 = []
 p = [:X₀ => 1.0, :p => 1.0, :d => 0.5]
 oprob = ODEProblem(rn, u0, tspan, p)
@@ -262,12 +263,12 @@ end
 ```
 Now, we can also declare our initial conditions and parameter values as vectors as well:
 ```@example dsl_advanced_vector_variables
-using OrdinaryDiffEqTsit5, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [:X => [0.0, 2.0]]
 tspan = (0.0, 1.0)
 ps = [:k => [1.0, 2.0]]
 oprob = ODEProblem(two_state_model, u0, tspan, ps)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 
@@ -489,12 +490,12 @@ X
 ```
 Next, we can now use these to e.g. designate initial conditions and parameter values for model simulations:
 ```@example dsl_advanced_programmatic_unpack
-using OrdinaryDiffEqTsit5, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [X => 0.1]
 tspan = (0.0, 10.0)
 ps = [p => 1.0, d => 0.2]
 oprob = ODEProblem(bd_model, u0, tspan, ps)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 plot(sol)
 ```
 

docs/src/model_creation/examples/basic_CRN_library.md

@@ -18,9 +18,9 @@ nothing # hide
 ```
 We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:
 ```@example crn_library_birth_death
-using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqDefault
 oprob = ODEProblem(bd_process, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 nothing # hide
 ```
 Next, a chemical Langevin equation-based SDE simulation:
@@ -52,7 +52,7 @@ plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))
 ## [Two-state model](@id basic_CRN_library_two_states)
 The two-state model describes a component (here called $X$) which can exist in two different forms (here called $X₁$ and $X₂$). It switches between these forms at linear rates. First, we simulate the model using both ODEs and SDEs:
 ```@example crn_library_two_states
-using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, StochasticDiffEq, Plots
 two_state_model = @reaction_network begin
     (k₁,k₂), X₁ <--> X₂
 end
@@ -62,7 +62,7 @@ tspan = (0.0, 1.0)
 ps = [:k₁ => 2.0, :k₂ => 3.0]
 
 oprob = ODEProblem(two_state_model, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
 
 sprob = SDEProblem(two_state_model, u0, tspan, ps)
@@ -96,9 +96,9 @@ u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
 tspan = (0., 100.)
 ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
 
-using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqDefault
 oprob = ODEProblem(mm_system, u0, tspan, ps)
-osol  = solve(oprob, Tsit5())
+osol  = solve(oprob)
 
 using StochasticDiffEq
 sprob = SDEProblem(mm_system, u0, tspan, ps)
@@ -132,14 +132,14 @@ end
 ```
 First, we perform a deterministic ODE simulation:
 ```@example crn_library_sir
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:S => 99, :I => 1, :R => 0]
 tspan = (0.0, 500.0)
 ps = [:α => 0.001, :β => 0.01]
 
 # Solve ODEs.
 oprob = ODEProblem(sir_model, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))
 ```
 Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
@@ -179,14 +179,14 @@ Below, we perform a simple deterministic ODE simulation of the system. Next, we
 
 In two separate plots.
 ```@example crn_library_cc
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
 tspan = (0., 15.)
 ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
 
 # solve ODEs
 oprob = ODEProblem(cc_system, u0, tspan, ps)
-osol  = solve(oprob, Tsit5())
+osol  = solve(oprob)
 
 plt1 = plot(osol; idxs = [:S₁, :S₂, :P], title = "Substrate and product dynamics")
 plt2 = plot(osol; idxs = [:C, :S₁C, :CP], title = "Catalyst and intermediaries dynamics")
@@ -206,16 +206,16 @@ end
 ```
 We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.
 ```@example crn_library_wilhelm
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0_1 = [:X => 1.5, :Y => 0.5]
 u0_2 = [:X => 2.5, :Y => 0.5]
 tspan = (0., 10.)
 ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
 
 oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
 oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
-osol1 = solve(oprob1, Tsit5())
-osol2 = solve(oprob2, Tsit5())
+osol1 = solve(oprob1)
+osol2 = solve(oprob2)
 plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
 plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
@@ -234,13 +234,13 @@ A simple example of such a loop is a transcription factor which activates its ow
 
 We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
 ```@example crn_library_self_activation
-using JumpProcesses, OrdinaryDiffEqTsit5, Plots
+using JumpProcesses, OrdinaryDiffEqDefault, Plots
 u0 = [:X => 4]
 tspan = (0.0, 1000.0)
 ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
 
 oprob = ODEProblem(sa_loop, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 
 jinput = JumpInputs(sa_loop, u0, tspan, ps)
 jprob = JumpProblem(jinput)
@@ -267,16 +267,16 @@ end
 ```
 It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst [automatically perform these adjustments](@ref introduction_to_catalyst_ratelaws), and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
 ```@example crn_library_brusselator
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X => 1.0, :Y => 1.0]
 tspan = (0., 50.)
 ps1 = [:A => 1.0, :B => 1.0]
 ps2 = [:A => 1.0, :B => 1.8]
 
 oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
 oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
-osol1  = solve(oprob1, Tsit5())
-osol2  = solve(oprob2, Tsit5())
+osol1  = solve(oprob1)
+osol2  = solve(oprob2)
 oplt1 = plot(osol1; title = "No Oscillation")
 oplt2 = plot(osol2; title = "Oscillation")
 
@@ -296,16 +296,16 @@ end
 ```
 Whether the Repressilator oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
 ```@example crn_library_brusselator
-using OrdinaryDiffEqTsit5, StochasticDiffEq, Plots
+using OrdinaryDiffEqDefault, StochasticDiffEq, Plots
 u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
 tspan = (0., 200.)
 ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
 ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]
 
 oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
 oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
-osol1  = solve(oprob1, Tsit5())
-osol2  = solve(oprob2, Tsit5())
+osol1  = solve(oprob1)
+osol2  = solve(oprob2)
 oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
 oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")
 
@@ -337,12 +337,12 @@ end
 ```
 Here we simulate the model for a single initial condition, showing both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
 ```@example crn_library_chaos
-using OrdinaryDiffEqTsit5, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
 tspan = (0.0, 50.0)
 p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7]
 oprob = ODEProblem(wr_model, u0, tspan, p)
-sol = solve(oprob, Tsit5())
+sol = solve(oprob)
 
 plt1 = plot(sol; title = "Time-state space")
 plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")

docs/src/model_creation/examples/hodgkin_huxley_equation.md

@@ -18,7 +18,7 @@ complete model.
 
 We begin by importing some necessary packages:
 ```@example hh1
-using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEqRosenbrock
+using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEqDefault
 ```
 
 ## Building the model via the Catalyst DSL
@@ -127,7 +127,7 @@ amplitude of the stimulus is non-zero and see if we get action potentials
 tspan = (0.0, 50.0)
 @unpack V,I₀ = hhmodel
 oprob = ODEProblem(hhmodel, u_ss, tspan, [I₀ => 10.0])
-sol = solve(oprob, Rosenbrock23())
+sol = solve(oprob)
 plot(sol, idxs = V, legend = :outerright)
 ```
 
@@ -188,7 +188,7 @@ eliminate the algebraic equations for the ionic currents when constructing the
 ```@example hh1
 @unpack I₀,V = hhmodel2
 oprob = ODEProblem(hhmodel2, u_ss, tspan, [I₀ => 10.0]; structural_simplify = true)
-sol = solve(oprob, Rosenbrock23())
+sol = solve(oprob)
 plot(sol, idxs = V, legend = :outerright)
 ```