BioPreDynB4.jmd

---
title: BioPreDyn-bench B4 Work-Precision Diagrams
author: Chris Rackauckas
---

The following benchmark is of 35 ODEs modeling the metabolism of Chinese
Hamster Ovary (CHO) cells from the
[BioPreDyn-bench](https://doi.org/10.1186/s12918-015-0144-4) benchmark suite
(Problem B4). The model describes a batch fermentation process with 32
reactions including glycolysis, TCA cycle, amino acid metabolism, and product
protein formation. The kinetics use log-linear elasticity coefficients with
117 parameters.

The model is manually translated from the original BioPreDyn-bench Matlab/C
files since the SBML representation uses rate rules rather than reactions.

```julia
using OrdinaryDiffEq, DiffEqDevTools, Sundials, Plots,
      LSODA, LinearAlgebra, BenchmarkTools,
      SciMLBenchmarks

gr()
```

## Load the model

```julia
include(joinpath(@__DIR__, "Models", "BioPreDyn", "b4_model.jl"))

prob = ODEProblem{true, SciMLBase.FullSpecialize}(
    biopredyn_b4!, B4_U0, B4_TSPAN, B4_NOMINAL_PARAMETERS)
@show length(prob.u0)
```

## Picture of the solution

```julia
sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
plot(sol; legend = false, fmt = :png,
     title = "BioPreDyn B4: CHO Cell Metabolism (35 ODEs)")
```

## Generate Reference Solution

```julia
@time sol = solve(prob, Rodas5P(), abstol = 1/10^14, reltol = 1/10^14)
test_sol = TestSolution(sol);
```

## Setups

```julia
default(legendfontsize = 7, framestyle = :box, gridalpha = 0.3, gridlinewidth = 2.5)
```

## Work-Precision Diagrams (High Tolerances)

```julia
abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4)

setups = [
    Dict(:alg => Rosenbrock23()),
    Dict(:alg => TRBDF2()),
    Dict(:alg => FBDF()),
    Dict(:alg => QNDF()),
    Dict(:alg => Rodas4()),
    Dict(:alg => Rodas5()),
    Dict(:alg => RadauIIA5()),
    Dict(:alg => lsoda()),
    Dict(:alg => CVODE_BDF()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
    saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["Rosenbrock23" "TRBDF2" "FBDF" "QNDF" "Rodas4" "Rodas5" "RadauIIA5" "lsoda" "CVODE_BDF"]
plot(wp; label = names, title = "BioPreDyn B4: High Tolerances")
```

## Work-Precision Diagrams (Low Tolerances)

```julia
abstols = 1.0 ./ 10.0 .^ (7:12)
reltols = 1.0 ./ 10.0 .^ (4:9)

setups = [
    Dict(:alg => FBDF()),
    Dict(:alg => QNDF()),
    Dict(:alg => Rodas4()),
    Dict(:alg => Rodas5P()),
    Dict(:alg => KenCarp4()),
    Dict(:alg => KenCarp47()),
    Dict(:alg => Kvaerno5()),
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => lsoda()),
    Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
    saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["FBDF" "QNDF" "Rodas4" "Rodas5P" "KenCarp4" "KenCarp47" "Kvaerno5" "CVODE_BDF" "lsoda" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: Low Tolerances")
```

## Work-Precision Diagrams (SDIRK Methods)

```julia
abstols = 1.0 ./ 10.0 .^ (5:10)
reltols = 1.0 ./ 10.0 .^ (2:7)

setups = [
    Dict(:alg => KenCarp3()),
    Dict(:alg => KenCarp4()),
    Dict(:alg => KenCarp5()),
    Dict(:alg => Kvaerno3()),
    Dict(:alg => Kvaerno4()),
    Dict(:alg => Kvaerno5()),
    Dict(:alg => TRBDF2()),
    Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
    saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["KenCarp3" "KenCarp4" "KenCarp5" "Kvaerno3" "Kvaerno4" "Kvaerno5" "TRBDF2" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: SDIRK Methods")
```

## Summary of Results

```julia
abstols = 1.0 ./ 10.0 .^ (5:12)
reltols = 1.0 ./ 10.0 .^ (2:9)

setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => QNDF()),
    Dict(:alg => FBDF()),
    Dict(:alg => Rodas4()),
    Dict(:alg => Rodas5P()),
    Dict(:alg => KenCarp4()),
    Dict(:alg => lsoda()),
    Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
    saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["CVODE_BDF" "QNDF" "FBDF" "Rodas4" "Rodas5P" "KenCarp4" "lsoda" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: Summary",
    left_margin = 10Plots.mm, right_margin = 10Plots.mm,
    legendfontsize = 10, tickfontsize = 10, guidefontsize = 10,
    legend = :topright, lw = 2, markersize = 6, size = (800, 600))
```

```julia, echo = false
using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
```