Systems

gene_expression_system(; mean_s=50, mean_x=mean_s, corr_time_s=1.0, corr_time_x=0.1, u0=SA[mean_s, mean_x], dtimes=0:0.1:2.0)

Creates a system for a very simple model of gene expression.

Model Description

This model consists of four reactions:

                κ
reaction 1:  ∅ ---> S

                λ
reaction 2:  S ---> ∅

                ρ
reaction 3:  S ---> S + X

                μ
reaction 4:  X ---> ∅

The first two reactions specify the evolution of the input signal S, and last two reactions specify the evolution of the output X. Thus, both the input and output signal are modeled as a simple birth-death process, however the birth rate of X increases with higher copy numbers of S.

Examples

The values of the reaction rates can be specified directly as follows:

using PathWeightSampling
system = PathWeightSampling.gene_expression_system(kappa = 10.0, lambda = 0.1, rho = 1.0, mu = 1.0)

# output

SimpleSystem with 4 reactions
Input variables: S(t)
Output variables: X(t)
Initial condition:
    S(t) = 50
    X(t) = 50
Parameters:
    κ = 10.0
    λ = 0.1
    ρ = 1.0
    μ = 1.0

Alternatively, the reaction rates can be specified indirectly through the following arguments:

• mean_s: The average copy number of S when the system has relaxed to steady state.
• mean_x: The average copy number of X when the system has relaxed to steady state.
• corr_time_s: The input signal correlation time. The shorter this time is, the faster the fluctuations in the input signal.
• corr_time_x: The output signal correlation time. This sets the timescale of output fluctuations.

using PathWeightSampling
system = PathWeightSampling.gene_expression_system(mean_s=25, mean_x=50, corr_time_s=1.0, corr_time_x=0.75)

# output

SimpleSystem with 4 reactions
Input variables: S(t)
Output variables: X(t)
Initial condition:
    S(t) = 25
    X(t) = 50
Parameters:
    κ = 25.0
    λ = 1.0
    ρ = 2.666666666666666
    μ = 1.3333333333333333

cooperative_chemotaxis_system(;
    lmax = 3,
    mmax = 9,
    n_clusters = 100,
    n_chey = 10_000,

    mean_l = 50,
    tau_l = 1.0,

    phosphorylate = 2000.0 / (n_chey * n_clusters),
    dephosphorylate = 2000.0 / (n_chey),
    dtimes = 0:0.1:20.0,
    varargs...
)

Create a system for a complex chemotaxis model.

Model Description

This model describes the bacterial chemotaxis signaling network.

Rml -> Rm(l+1), with rate kon L(t) (lmax - l): ligand binding to active state
Rm(l+1) -> Rml, with rate koff_A l: ligand unbinding
Rml -> R(m+1)l with rate km: methylation rate
R(m+1)l -> Rml with rate kdm: demethylation rate

Examples

Create a chemotaxis system with default parameters.

using PathWeightSampling
PathWeightSampling.cooperative_chemotaxis_system()

# output

ComplexSystem with 175 reactions
Input variables: L(t)
Latent variables: R_0_0(t), R_0_1(t), R_0_2(t), R_0_3(t), R_0_4(t), R_0_5(t), R_0_6(t), R_0_7(t), R_0_8(t), R_0_9(t), R_1_0(t), R_1_1(t), R_1_2(t), R_1_3(t), R_1_4(t), R_1_5(t), R_1_6(t), R_1_7(t), R_1_8(t), R_1_9(t), R_2_0(t), R_2_1(t), R_2_2(t), R_2_3(t), R_2_4(t), R_2_5(t), R_2_6(t), R_2_7(t), R_2_8(t), R_2_9(t), R_3_0(t), R_3_1(t), R_3_2(t), R_3_3(t), R_3_4(t), R_3_5(t), R_3_6(t), R_3_7(t), R_3_8(t), R_3_9(t)
Output variables: Yp(t), Y(t)
Initial condition:
    L(t) = 50
    R_0_0(t) = 100
    R_0_1(t) = 0
    R_0_2(t) = 0
    R_0_3(t) = 0
    R_0_4(t) = 0
    R_0_5(t) = 0
    R_0_6(t) = 0
    R_0_7(t) = 0
    R_0_8(t) = 0
    R_0_9(t) = 0
    R_1_0(t) = 0
    R_1_1(t) = 0
    R_1_2(t) = 0
    R_1_3(t) = 0
    R_1_4(t) = 0
    R_1_5(t) = 0
    R_1_6(t) = 0
    R_1_7(t) = 0
    R_1_8(t) = 0
    R_1_9(t) = 0
    R_2_0(t) = 0
    R_2_1(t) = 0
    R_2_2(t) = 0
    R_2_3(t) = 0
    R_2_4(t) = 0
    R_2_5(t) = 0
    R_2_6(t) = 0
    R_2_7(t) = 0
    R_2_8(t) = 0
    R_2_9(t) = 0
    R_3_0(t) = 0
    R_3_1(t) = 0
    R_3_2(t) = 0
    R_3_3(t) = 0
    R_3_4(t) = 0
    R_3_5(t) = 0
    R_3_6(t) = 0
    R_3_7(t) = 0
    R_3_8(t) = 0
    R_3_9(t) = 0
    Yp(t) = 0
    Y(t) = 10000
Parameters:
    κ = 50.0
    λ = 1.0
    E0 = 3.0
    δg = 2.995732273553991
    δf = -1.5
    lba = 0.05
    lbi = 0.05
    lda = 25.0
    ldi = 1.25
    mda = 0.03333333333333333
    mbi = 0.03333333333333333
    μ = 0.2
    ρ = 0.002

SimpleSystem(input_network, output_network, u0, ps, px, dtimes)

A SimpleSystem is a system used to compute the mutual information when there are no latent variables that need to be integrated out.

A SimpleSystem consists of two of ModelingToolkit's ReactionSystems: One that generates the input trajectories (input_network), and one that generates the output trajectories (output_network).

u0 represents the initial condition of the system.

ps and px are the parameter vectors of the input and output systems, respectively.

Examples

Create a SimpleSystem from a set of coupled birth-death processes.

using PathWeightSampling, Catalyst, StaticArrays

sn = @reaction_network begin
    κ, ∅ --> 2L
    λ, L --> ∅
end κ λ

xn = @reaction_network begin
    ρ, L --> L + 2X
    μ, X --> ∅
end ρ μ

u0 = SA[10, 20]
dtimes = 0:0.5:10.0
ps = [5.0, 1.0]
px = [3.0, 0.1]

system = PathWeightSampling.SimpleSystem(sn, xn, u0, ps, px, dtimes)

# output

SimpleSystem with 4 reactions
Input variables: L(t)
Output variables: X(t)
Initial condition:
    L(t) = 10
    X(t) = 20
Parameters:
    κ = 5.0
    λ = 1.0
    ρ = 3.0
    μ = 0.1

ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes, dist=nothing; aggregator=Direct())

A ComplexSystem is a system used to compute the mutual information with latent variables that need to be integrated out.

Therefore, a ComplexSystem consists of three ModelingToolkit ReactionSystems: One that generates the input trajectories (sn), one that models the latent variables (rn), and one that generates the output trajectories (xn).

Examples

using PathWeightSampling, Catalyst, StaticArrays

sn = @reaction_network begin
    κ, ∅ --> 2L
    λ, L --> ∅
end κ λ

rn = @reaction_network begin
    ρ, L + R --> L + LR
    μ, LR --> R
    ξ, R + CheY --> R + CheYp
    ν, CheYp --> CheY
end ρ μ ξ ν

xn = @reaction_network begin
    δ, CheYp --> CheYp + X
    χ, X --> ∅
end δ χ

u0 = SA[10, 30, 0, 50, 0, 0]
dtimes = 0:0.5:10.0
ps = [5.0, 1.0]
pr = [1.0, 4.0, 1.0, 2.0]
px = [1.0, 1.0]

system = PathWeightSampling.ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes)

# output

ComplexSystem with 8 reactions
Input variables: L(t)
Latent variables: R(t), LR(t), CheY(t), CheYp(t)
Output variables: X(t)
Initial condition:
    L(t) = 10
    R(t) = 30
    LR(t) = 0
    CheY(t) = 50
    CheYp(t) = 0
    X(t) = 0
Parameters:
    κ = 5.0
    λ = 1.0
    ρ = 1.0
    μ = 4.0
    ξ = 1.0
    ν = 2.0
    δ = 1.0
    χ = 1.0