Make more clear that mRNA has a synthesis term

[ftavella]

In the current GRN implementation we have:

# mRNA
push!(rxs, Reaction(1 - α[i]*m[i], nothing, [m[i]]))

which obscures the fact that mRNA has both a synthesis and degradation term. We could make this more clear by separating it into two lines:

# mRNA
push!(rxs, Reaction(1, nothing, [m[i]]))
push!(rxs, Reaction(α[i],  [m[i]], nothing))

[biphy]

The following is one possible demonstration of their identity based on the SciML API. It shows the two implementations are the same in terms of the symbolic ODEs they generate.

using Catalyst
using ModelingToolkit

using BiologicalOscillations

function modified_gene_regulatory_network(connectivity::AbstractMatrix)
    # Check that input is correct
    is_valid, errmsg = is_valid_connectivity(connectivity)
    if !is_valid
        throw(DomainError(connectivity, errmsg))
    end
    # Number of nodes
    N = size(connectivity, 1) 
    # Number of edges
    E = count(!iszero, connectivity)
    @variables t
    @species (m(t))[collect(1:N)]
    @species (X(t))[collect(1:N)]
    @parameters α[1:N], β[1:N], δ[1:N], γ[1:E], κ[1:E], η[1:E]
    rxs = Reaction[]
    # Node-specific reactions for mRNA and Protein
    for i=1:size(connectivity, 1)
        # mRNA

        # >>> modification here >>>
        push!(rxs, Reaction(1, nothing, [m[i]]))
        push!(rxs, Reaction(α[i], [m[i]], nothing))
        # <<< modification here <<<

        # Protein
        push!(rxs, Reaction(β[i]*m[i], nothing, [X[i]]))
        push!(rxs, Reaction(δ[i], [X[i]], nothing))
    end
    # Reactions between nodes
    e = 1 
    for (i, row) in enumerate(eachrow(connectivity))
        for (j, val) in enumerate(row)
        if val == -1
            # Repression
            push!(rxs, Reaction(hillr(abs(X[j]),γ[e],κ[e],η[e]), nothing, [m[i]]))
            e += 1
        elseif val == 1
            # Activation
            push!(rxs, Reaction(hill(abs(X[j]),γ[e],κ[e],η[e]), nothing, [m[i]]))
            e += 1
        end
        end
    end
    @named model = ReactionSystem(rxs, t)
    return model
end

connectivity1 = [1 0 -1 ; 1 0 0 ; 0 1 0]
connectivity2 = [0 0 0 -1 ; -1 1 0 0 ; 0 1 0 -1 ; 0 1 -1 0]

# system 1, original vs modified
rsys1 = BiologicalOscillations.gene_regulatory_network(connectivity1)
rsys_modified1 = modified_gene_regulatory_network(connectivity1)

# system 2, original vs modified
rsys2 = BiologicalOscillations.gene_regulatory_network(connectivity2)
rsys_modified2 = modified_gene_regulatory_network(connectivity2)

# convert defined in Catalyst.jl
op1 = convert(ODESystem, rsys1)
op_modified1 = convert(ODESystem, rsys_modified1)
op2 = convert(ODESystem, rsys2)
op_modified2 = convert(ODESystem, rsys_modified2)

# For the same system, the original and modified implementations must be equivalent
@assert ModelingToolkit.get_eqs(op1) == ModelingToolkit.get_eqs(op_modified1)
@assert ModelingToolkit.get_eqs(op2) == ModelingToolkit.get_eqs(op_modified2)

# Different systems
@assert ModelingToolkit.get_eqs(op1) != ModelingToolkit.get_eqs(op2)
@assert ModelingToolkit.get_eqs(op1) != ModelingToolkit.get_eqs(op_modified2)
@assert ModelingToolkit.get_eqs(op_modified1) != ModelingToolkit.get_eqs(op2)
@assert ModelingToolkit.get_eqs(op_modified1) != ModelingToolkit.get_eqs(op_modified2)

[ftavella]

This looks great to me!