Implement Matlab/Simbiology concepts (e.g. compartments, repeated assignments) in Julia/Catalyst

valerdam · May 14, 2024, 6:30pm

I am testing ways to code Matlab/Simbiology concepts (e.g. compartments, repeated assignments) in Julia/Catalyst and I would appreciate both general pointers on where to find more information as well as figuring out one piece of this puzzle in the code below. (error pasted as comment)

# test example to create a simple model with two compartments

using Catalyst
using DifferentialEquations 
using Plots   

@variables t

@parameters p = 1

# define compartment c1 having the volume c1
@parameters c1 = 10     
c1rn = @reaction_network c1 

# define compartment c1 having the volume c1
@parameters c2 = 2
c2rn = @reaction_network c2 

# define species D in compartment c1
@species    c1₊D(t) = 300
# define species D in compartment c2 
@species    c2₊D(t) = 100

# define species Dc in compartments c1 and c2 to represent concentrations
@variables  c1₊Dc(t)
@variables  c2₊Dc(t)

# define the reactions  

# This does not works as it uses the c1₊Dc and c2₊Dc variables in the rates 
rxn = [
    @reaction p*c1₊Dc, c1₊D => c2₊D
    @reaction p*c2₊Dc, c2₊D => c1₊D
]

#= This works but it defeats the purpose of creating variables to make the definition of the reaction more readable
rxn = [
    @reaction p*c1₊D/c1, c1₊D => c2₊D
    @reaction p*c2₊D/c2, c2₊D => c1₊D
]
=#

# define the equations defining c1₊Dc and c2₊Dc used in the rates above
ra = [ 
	c1₊Dc ~ c1₊D / c1
	c2₊Dc ~ c2₊D / c2
]

# build the reaction system with component systems c1rn and c2rn corresponding to the two compartments
@named rs = ReactionSystem([rxn; ra], t, [c1₊D, c2₊D, c1₊Dc, c2₊Dc], [p, c1, c2]; systems=[c1rn,c2rn])

meq0 = convert(ODESystem, rs)
meq = structural_simplify(meq0)
tspan = (0.0,250.0)

ode = ODEProblem(meq, [], tspan, [])
# all works up to this point
# but fails on the following line with 
# LoadError: UndefVarError: `c1₊Dc` not defined

sol = solve(ode, Tsit5())
plot(sol)

isaacsas · May 15, 2024, 1:21am

What about this:

using Catalyst, DifferentialEquations, Plots

t = Catalyst.DEFAULT_IV

c1 = @reaction_network c1 begin
    @parameters c1 = 10
    @species D(t) = 300
    @variables Dc(t)
end

c2 = @reaction_network c2 begin
    @parameters c2 = 10
    @species D(t) = 300
    @variables Dc(t)
end

@parameters p = 1
rxs =  [Reaction(p*c1.Dc, [c1.D], [c2.D]),
        Reaction(p*c2.Dc, [c2.D], [c1.D]),
        (c1.Dc ~ c1.D / c1.c1),
        (c2.Dc ~ c2.D / c2.c2)]
@named rs = ReactionSystem(rxs, t; systems = [c1, c2])

osys = convert(ODESystem, rs)
osys = structural_simplify(osys)

tspan = (0.0,250.0)
ode = ODEProblem(osys, [], tspan, [])
sol = solve(ode, Tsit5())
plot(sol)

Note though, the way you are specifying reactions you are generating a nonlinear reaction, i.e. your final ODEs are:

julia> full_equations(osys)
2-element Vector{Equation}:
 Differential(t)(c1₊D(t)) ~ (-p*(c1₊D(t)^2)) / c1₊c1 + (p*(c2₊D(t)^2)) / c2₊c2
 Differential(t)(c2₊D(t)) ~ (p*(c1₊D(t)^2)) / c1₊c1 + (-p*(c2₊D(t)^2)) / c2₊c2

This is because Catalyst automatically includes the mass action portion of the rate from the substrates of the reaction and multiplies the rest of the rate by it. So perhaps you want something more like

using Catalyst, DifferentialEquations, Plots

t = Catalyst.DEFAULT_IV

c1 = @reaction_network c1 begin
    @parameters c1 = 10
    @species D(t) = 300
end

c2 = @reaction_network c2 begin
    @parameters c2 = 10
    @species D(t) = 300
end

@parameters p = 1
rxs =  [Reaction(p/c1.c1, [c1.D], [c2.D]),
        Reaction(p/c2.c2, [c2.D], [c1.D])]
@named rs = ReactionSystem(rxs, t; systems = [c1, c2])

osys = convert(ODESystem, rs)
osys = structural_simplify(osys)
tspan = (0.0,250.0)
ode = ODEProblem(osys, [], tspan, [])
sol = solve(ode, Tsit5())
plot(sol)

so that

julia> full_equations(osys)
2-element Vector{Equation}:
 Differential(t)(c1₊D(t)) ~ (-p*c1₊D(t)) / c1₊c1 + (p*c2₊D(t)) / c2₊c2
 Differential(t)(c2₊D(t)) ~ (p*c1₊D(t)) / c1₊c1 + (-p*c2₊D(t)) / c2₊c2

If you prefer to specify the full rate yourself and not have Catalyst infer the portion arising from the substrates you can use the only_use_rate = true keyword argument, see here. This is generally not recommended however, for example resulting in less efficient Gillespie simulations if you want to look at stochastic models (less efficient because the reaction is no longer treated as a basic mass action reaction and instead uses a more general representation internally).

isaacsas · May 15, 2024, 1:36am

We’ve discussed adding compartments before, but haven’t settled on what a good design would be. Should a compartment store one associated ReactionSystem, should it be an optional component within a ReactionSystem, etc?

For now I would suggest building one ReactionSystem for each fixed compartment, and then building a bridge system to handling coupling them together (i.e. for transport between compartments).

isaacsas · May 15, 2024, 1:38am

For references, probably the compositional modeling tutorial would be useful here.

Note that we are currently revamping the docs and making a number of related updates to the DSL to simplify model creation that would make this example simpler I believe. These should be released in the near future as part of Catalyst 14.

valerdam · May 15, 2024, 3:05am

Many thanks for your prompt response.

While it gave me some great pointers it does not really solve the issue. Apologies - I may not have articulated my question very clearly.

I do not believe your point about nonlinear reactions is correct. In the definition of the reactions I have used => instead of → which is supposed to instruct Catalyst to not generate mass action and use the rate as I give. In my version uncommenting out the portion of the code that labeled as working the full equations generated are linear. I did not realize you can achieve the same with only_use_rate = true - thanks for the tip.
I do want my reactions to be more general than mass action. For the systems I am interested in solving the rates are quite massive nonlinear expressions.
My main problem was in defining the intermediate variables c1₊Dc. In the example here they are trivial in my use case they can be massive expressions that would make the code impossible to maintain. I see they work fine if they are used only as outputs but i my case I need to refer to them in in the rate expressions for other reactions. This is what I was trying to achieve with my code. It seems that while the substitutions are being kept in the model equations (meq in my code) but they are not used to make the appropriate substitutions in the equations which still make reference to the intermediate variables.
Full disclosure - I have been working with Matlab and Simbiology for decades but only started with Julia. In Simbiology there is a concept of “repeated assignment” which if I were to wright the RHS myself are just assignments executed in the right order and then used in the rate equations that follow - like local variables in the RHS. Some of these may also be outputs but that is already taken care of nicely.

Again many thanks for your reply.

valerdam · May 15, 2024, 3:41am

I am biased here due to my long experience using Matlab/Simbiology, but perhaps my comments below can be of help if you plan to implement a compartment concept in Catalyst.

In Simbiology compartments have a volume and contain species that can be expressed either as amount or concentration. Volumes can also be modified by equations or forcing functions (e.g growth of the volume of a cell).
Reactions belong to the model not to the compartment. There are quite a number of rate expressions for species transport from one compartment to the other. Simbiology takes care of scaling the rate by the compartment volumes as appropriate. It is also possible for a species in one compartment to be used in an equation in another (e.g a receptor siting on the cell membrane can influence reactions happening in the cell or in the medium, or influencing the transport of other species).
Parameters can be local to a compartment, local to a reaction or global to the model. Parameters can also be constant or variable - i.e. they can be defined by the equations.
There are several types of rules - most used types are: Initial assignments - which you are already implemented nicely with the defaults. Repeated assignments - as I mentioned in my example they are effectively intermediate variables used in constructing the rete equations. Rate rules - which correspond to defining explicitly differential equations - this is already available.
As far as I know there isn’t yet an explicit submodel concept (as in Simulink). The compositional model here seem to be more alike the Simulink one.
There are also units. Every species/parameter/ compartment volume has a unit and if consistent but not commensurate the unit conversion is done automatically.
Finally there are doses. This are quantities that are added to specific species at predefined times. One can achieve these in Julia using events but in Simbiology they are first class citizens.
Please let me know if this is helpful.

isaacsas · May 15, 2024, 10:44am

I missed you used =>, so yes that will give the rate law you want.

My main problem was in defining the intermediate variables c1₊Dc. In the example here they are trivial in my use case they can be massive expressions that would make the code impossible to maintain. I see they work fine if they are used only as outputs but i my case I need to refer to them in in the rate expressions for other reactions. This is what I was trying to achieve with my code. It seems that while the substitutions are being kept in the model equations (meq in my code) but they are not used to make the appropriate substitutions in the equations which still make reference to the intermediate variables.

I’m not sure what you mean by this. Do you have another example of the problem you are mentioning now? Doesn’t my first example let you use the algebraic variables within the reactions as you want (modulo that you need to add only_use_rate = true in the reaction construction)?

Parameters can be local to a compartment, local to a reaction or global to the model. Parameters can also be constant or variable - i.e. they can be defined by the equations.

Parameters as functions of other parameters should work, i.e.

@parameters begin
    k1 
    k2 = 2*k1
end

They can’t be functions of species/variables (because at that point they should be a variable themselves).

There are also units. Every species/parameter/ compartment volume has a unit and if consistent but not commensurate the unit conversion is done automatically.

Units should be supported and checked, but conversions are not handled automatically. This is admittedly a feature that few have used, so may not cover everything one would want – feel free to open issues on the Catalyst Github if you find problems using units or are looking for missing functionality.

Finally there are doses. This are quantities that are added to specific species at predefined times. One can achieve these in Julia using events but in Simbiology they are first class citizens.

Do you mean that you’d like to have “doses” as an explicit model component instead of having to add symbolic events to model the doses within a ReactionSystem?

isaacsas · May 15, 2024, 11:10am

Just for reference, this is how I’d setup your example with Catalyst 14 (which isn’t released yet but will hopefully be done this month):

using Catalyst, DifferentialEquations, Plots

t = default_t()

c1 = @network_component c1 begin
    @parameters V = 10
    @species D(t) = 300
    @variables Dc(t)
    @equations Dc ~ D / V
end

c2 = @network_component c2 begin
    @parameters V = 10
    @species D(t) = 300
    @variables Dc(t)
    @equations Dc ~ D / V
end

@parameters p = 1
rxs =  [Reaction(p/c1.V, [c1.D], [c2.D]; use_only_rate = true),
        Reaction(p/c2.V, [c2.D], [c1.D]; use_only_rate = true)]
@named rs = ReactionSystem(rxs, t; systems = [c1, c2])
rs = complete(rs)

osys = convert(ODESystem, rs)
osys = structural_simplify(osys)

tspan = (0.0,250.0)
ode = ODEProblem(osys, [], tspan, [])
sol = solve(ode, Tsit5())
plot(sol)

valerdam · May 15, 2024, 1:58pm

Many thanks for all your comments.

The very first version adding ‘only_use_rate = true’ works as I intended. (when I tried it first I did not modify it and it did not gave me the results I was expecting).

The second version removes the c1.Dc and c2.Dc variables and that was the main reason for doing this test. Please disregard my comment about “intermediate variables” as your first version of the code addresses the issue.

On the doses topic: Doses are widely used in pharmacology domain applications. They can be implemented using events but it would be nice to encapsulate them in an explicit model component. Repeat doses (e.g daily doses) are common or infusion doses where the species is infused linearly.

isaacsas · May 15, 2024, 2:13pm

Probably the best way to handle dosing is to write a couple functions that generate events for the dose types you want, and then you can just call them as part of your code. MTK events do still need some work to be more flexible so that one can give a symbolic variable for the time of a discrete event, or specify directionality for when a continuous event can occur.

I’m not sure dosing makes sense to explicitly add to Catalyst, but perhaps such functionality would make sense as part of a higher-level library of useful Catalyst add-on components (i.e. standard network motifs, standard biological event types, etc – something like ModelingToolkitStandardLibrary but for Catalyst).