ReactiveDynamics.jl

About | Context | Four Sketches | Documentation

The package provides a category of reaction (transportation) network-type problems formalized on top of the generalized algebraic theory, and is compatible with the SciML ecosystem.

Our motivation stems from the area of system dynamics, which is a mathematical modeling approach to frame, analyze, and optimize complex (nonlinear) dynamical systems, to augment the strategy and policy design.

The central concept is of a transition (transport, flow, rule, reaction - a subcategory of general algebraic action). Generally, a transition prescribes a stochastic rule which repeatedly transforms the modeled system's resources. An elementary instance of such an ontology is provided by chemical reaction networks.

A reaction network (system modeled) is then a tuple $(T, R)$, where $T$ is a set of the transitions and $R$ is a set of the network's resource classes (aka species). The simultaneous action of transitions on the resources evolves the dynamical system.

The transitions are generally stateful (i.e., act over a period of time). Moreover, at each time step a quantum of the a transition's instances is brought into the scope, where the size of the batch is drived by a Poisson counting process. A transition takes the from rate, a*A + b*B + ... --> c*C + ..., where rate gives the expected batch size per time unit. A, B, etc., are the resources, and a, b, etc., are the generalized stoichiometry coefficients. Note that both rate and the "coefficients" can in fact be given by a function which depends on the system's instantaneous state (stochastic, in general). In particular, even the structural form of a transition can be stochastic, as will be demonstrated shortly.

An instance of a stateful transition evolves gradually from its genesis up to the terminal point, where the products on the instance's right hand-side are put into the system. An instance advances proportionally to the quantity of resources allocated. To understand this behavior, we note that the system's resource classes (or occurences of a resource class in a transition) may be assigned a modality; a modality governs the interaction between the resource and the transition, as well as the interpretation of the generalized stoichiometry coefficient.

In particular, a resource can be allocated to an instance either for the instance's lifetime or a single time step of the model's evolution (after each time step, the resource may be reallocated based on the global demand). Moreover, the coefficient of a resource on the left hand-side can either be interpreted as a total amount of the resource required or as an amount required per time unit. Similarly, it is possible to declare a resource conserved, in which case it is returned into the scope once the instance terminates.

The transitions are parametric. That is, it is possible to set the period over which an instance of a transition acts in the system (as well as the maximal period of this action), the total number of transition's instances allowed to exist in the system, etc. An annotated transition takes the form rate, a*A + b*B + ... --> c*C + ..., prm => val, ..., where the numerical values can be given by a function which depends on the system's state. Internally, the reaction network is represented as an attributed C-set.

For an overview of accepted attributes for both transitions and species classes, read the docs

A network's dynamics is specified using a compact modeling metalanguage. Moreover, we have integrated another expression comprehension metalanguage which makes it easy to generate arbitrarily complex dynamics from a single template transition!

Taking unions of reaction networks is fully supported, and it is possible to identify the resource classes as appropriate.

Moreover, it is possible to export and import reaction network dynamics using the TOML format.

Once a network's dynamics is specified, it can be converted to a problem and simulated. The exported problem is a DiscreteProblem compatible with DifferentialEquations.jl ecosystem, and hence the latter package's all powerful capabilities are available. For better user experience, we have tailored and exported many of the functionalities within the modeling metalanguage, including ensemble analysis, parameter optimization, parameter inference, etc. Would you guess that universal differential equations are supported? If even the dynamics is unknown, you may just infer it!

The package is an integral part of the Dynamics of Value Evolution (DyVE) computational framework for learning, designing, integrating, simulating, and optimizing R&D process models, to better inform strategic decisions in science and business.

As the framework evolves, multiple functionalities have matured enough to become standalone packages.

This includes AlgebraicAgents.jl, a lightweight package to enable hierarchical, heterogeneous dynamical systems co-integration. It implements a highly scalable, fully customizable interface featuring sums and compositions of dynamical systems. In present context, we note it can be used to co-integrate a reaction network problem with, e.g., a stochastic ordinary differential problem!

CEEDesigns.jl is a decision-making framework for the cost-efficient design of experiments, with direct applications in drug research, healthcare, and various other fields. Typically, a design consists of multiple experiments. Each experiment is regarded as an option to acquire additional experimental evidence and is associated with a monetary cost and an execution time. The framework, then, aims to select experiments that balance the value of acquired information and the incurred costs.

Another package is GeneratedExpressions.jl, a metalanguage to support code-less expression comprehensions. In the present context, expression comprehensions are used to generate complex dynamics from user-specified template transitions.

For other examples, see the tutorials.

The acronym SIR stands for susceptible, infected, and recovered, and as such the SIR model attempts to capture the dynamics of disease spread. We express the SIR dynamics as a reaction network using the compact modeling metalanguage.

Follow the SIR model's reactions:

using ReactiveDynamics

# model dynamics
sir_acs = @ReactionNetwork begin
        α*S*I, S+I --> 2I, name=>I2R
        β*I, I --> R, name=>R2S 
end

# simulation parameters
## initial values
@prob_init sir_acs S=999 I=10 R=0
## uncertainty in initial values (Gaussian)
@prob_uncertainty sir_acs S=10. I=5.
## parameters
@prob_params sir_acs α=0.0001 β=0.01
## other arguments passed to the solver
@prob_meta sir_acs tspan=250 dt=.1

The resulting reaction network is represented as an attributed C-set:

Next we solve the problem.

# turn model into a problem
prob = @problematize sir_acs

# solve the problem over multiple trajectories
sol = @solve prob trajectories=20

# plot the solution
@plot sol plot_type=summary
## show only species S
@plot sol plot_type=summary show=:S
## plot evolution over (0., 100.) in green (propagates to Plots.jl)
@plot sol plot_type=summary c=:green xlimits=(.0, 100.)

Before we move on to more intricate examples demonstrating generative capabilities of the package, let's sketch a toy pharma model with as little as three transitions.

toy_pharma_model = @ReactionNetwork

First, a "discovery" transition will take a team of scientist and a portion of a company's budget at the input (say, for experimental resources), and it will output candidate compounds.

@push toy_pharma_model α(candidate_compound, marketed_drug, κ) 3*@conserved(scientist) + @rate(budget) --> candidate_compound name=>discovery probability=>.3 cycletime=>6 priority=>.5

Note that per a time unit, α(candidate_compound, marketed_drug, κ) "discovery" projects will be started. We provide a name of the class of transitions (name=>discovery), set up a probability of the transition terminating successfully (probability=>.3), a cycle time (cycletime=>6), and we provide a weight of the transitions' class for use in resource allocation (priority=>.5).

Moreover, we annotate "scientists" as a conserved resource (no matter how the project terminates, the workforce isn't consumed), i.e., @conserved(scientist), and we state that a unit "budget" is consumed per a time unit, i.e., @rate(budget).

Next, candidate compounds will undergo clinical trials. If successful, a compound transforms into a marketed drug, and the company receives a premium.

@push toy_pharma_model β(candidate_compound, marketed_drug) candidate_compound + 5*@conserved(scientist) + 2*@rate(budget) --> marketed_drug + 5*budget name=>dx2market probability=>.5+.001*@t() cycletime=>4

Note that as time evolves, the probability of technical success increases, i.e., probability=>.5+.001*@t().

In addition, marketed drugs bring profit to the company - which will fuel invention of new drugs.

We model the situation as a periodic callback.

@periodic toy_pharma_model 1. budget += 11*marketed_drug

A marketed drug may eventually be withdrawn from the market. To account for such scenario, we add the following transition:

@push toy_pharma_model  γ*marketed_drug marketed_drug --> ∅ name=>drug_withdrawn

Next we provide the functions α and β.

@register α(number_candidate_compounds, number_marketed_drugs, κ) = κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
@register β(number_candidate_compounds, number_marketed_drugs) = numbercandidate_compounds + exp(-number_marketed_drugs)

Likewise, we set the remaining parameters, initial values, and model metadata:

# simulation parameters
## initial values
@prob_init toy_pharma_model candidate_compound=5 marketed_drug=6 scientist=20 budget=100
## parameters
@prob_params toy_pharma_model κ=4 γ=.1
## other arguments passed to the solver
@prob_meta toy_pharma_model tspan=250 dt=.1

And we problematize the model, solve the problem, and plot the solution:

prob = @problematize toy_pharma_model
sol = @solve prob trajectories=20
@plot sol plot_type=summary show=[:candidate_compound, :marketed_drug]

We introduce a complex reaction network as a union of $n_{\text{models}}$ reaction networks, where the off-diagonal interactions are sparse.

To harness the capabilities of GeneratedExpressions.jl, let us first declare a template atomic model.

# submodel.jl
# substitute $r as the global number of resources, $i as the submodel identifier
@register begin 
        push!(ns, rand(1:5)); ϵ = 10e-2
        push!(M, rand(ns[$i], ns[$i])); foreach(i -> M[$i][i, i] += ϵ, 1:ns[$i])
        foreach(i -> M[$i][i, :] /= sum(M[$i][i, :]), 1:ns[$i])
        push!(cycle_times, rand(1:5, ns[$i], ns[$i]))
        push!(demand, rand(1:10, ns[$i], ns[$i], $r)); push!(production, rand(1:10, ns[$i], ns[$i], $r))
end
    
# generate submodel dynamics
push!(rd_models, @ReactionNetwork begin
                M[$i][$m, $n], state[$m] + {demand[$i][$m, $n, $l]*resource[$l], l=1:$r, dlm=+} --> state[$n] + 
                        {production[$i][$m, $n, $l]*resource[$l], l=1:$r, dlm=+}, cycle_time=>cycle_times[$i][$m, $n], probability_of_success=>$m*$n/(n[$i])^2
        end m=1:ReactiveDynamics.ns[$i] n=1:ReactiveDynamics.ns[$i]
)

The next step is to instantiate the atomic models (submodels).

using ReactiveDynamics
## setup the environment
rd_models = ReactiveDynamics.ReactionNetwork[] # submodels

# needs to live within ReactiveDynamics's scope
# the arrays will contain the submodel
@register begin
    ns = Int[] # size of submodels
    M = Array[] # transition intensities of submodels
    cycle_times = Array[] # cycle times of transitions in submodels
    demand = Array[]; production = Array[] # resource production / generation for transitions in submodels
end

Load $n_{\text{models}}$ the atomic models (substituting into submodel.jl):

n_models = 5; r = 2 # number of submodels, resources

# submodels: dense interactions
@generate {@fileval(submodel.jl, i=$i, r=r), i=1:n_models}

We take union of the atomic models, and we identify common resources.

# batch join over the submodels
rd_model = @generate "@join {rd_models[\$i], i=1:n_models, dlm=' '}"

# identify resources
@generate {@equalize(rd_model, @alias(resource[$j])={rd_models[$i].resource[$j], i=1:n_models, dlm=:(=)}), j=1:r}

Next step is to add some off-diagonal interactions.

# sparse off-diagonal interactions, sparse declaration
# again, we use GeneratedExpressions.jl
sparse_off_diagonal = zeros(sum(ReactiveDynamics.ns), sum(ReactiveDynamics.ns))
for i in 1:n_models
    j = rand(setdiff(1:n_models, (i, )))
    i_ix = rand(1:ReactiveDynamics.ns[i]); j_ix = rand(1:ReactiveDynamics.ns[j])
    sparse_off_diagonal[i_ix+sum(ReactiveDynamics.ns[1:i-1]), j_ix+sum(ReactiveDynamics.ns[1:j-1])] += 1
    interaction_ex = """@push rd_model begin 1., var"rd_models[$i].state[$i_ix]" --> var"rd_models[$j]__state[$j_ix]" end"""
    eval(Meta.parseall(interaction_ex))
end

Let's plot the interactions:

The resulting model can then be conveniently simulated using the modeling metalanguage.

using ReactiveDynamics: nparts
u0 = rand(1:1000, nparts(rd_model, :S))
@prob_init rd_model u0

@prob_meta rd_model tspan=10

prob = @problematize rd_model
sol = @solve prob trajectories=2

# plot "state" species only
@plot sol plot_type=summary show=r"state"

We demonstrate how to fit unknown part of dynamics to empirical data.

We use @register to define a simple linear function within the scope of module ReactiveDynamics; parameters of the function will be then optimized for. Note that function_to_learn can be generally replaced with a neural network (Flux chain), etc.

## some embedded function (neural network, etc.)
@register begin
    function function_to_learn(A, B, C, params)
        [A, B, C]' * params # params: 3-element vector
    end
end

Next we set up a simple dynamics and supply initial parameters.

acs = @ReactionNetwork begin
    function_to_learn(A, B, C, params), A --> B+C
    1., B --> C
    2., C --> B
end

# initial values, check params
@prob_init acs A=60. B=10. C=150.
@prob_params acs params=[.01, .01, .01]
@prob_meta acs tspan=100.

Let's next see the numerical results for the initial guess.

sol = @solve acs
@plot sol

Next we supply empirical data and fit params.

time_points = [1, 50, 100]
data = [60 30 5]

@fit_and_plot acs data time_points vars=[A] params α maxeval=200 lower_bounds=0 upper_bounds=.01