Implement a split enzyme network model for topology exploration

.github/workflows/CI.yml

@@ -19,7 +19,6 @@ jobs:
       matrix:
         version:
           - '1.8'
-          - 'nightly'
         os:
           - ubuntu-latest
         arch:

docs/pages.jl

@@ -6,5 +6,5 @@ pages = Any["Home" => "index.md",
                                               "gene_regulatory_network/grn_model_details.md", 
                                               "gene_regulatory_network/grn_documentation.md"],
             "Feature Calculation" => Any["feature_calculation/fc_documentation.md"],
-            "Models" => Any["models/elowitz_2000.md", "models/guan_2008.md"],
+            "Models" => Any["models/elowitz_2000.md", "models/guan_2008.md", "models/kimchi_2020.md"],
 ]

docs/src/models/kimchi_2020.md

@@ -0,0 +1,34 @@
+# [Kimchi et al. 2020 - Posttranslational oscillator](@id kimchi_2020)
+```@docs
+kimchi_2020
+```
+```@setup kimchi_2020
+using BiologicalOscillations, Catalyst, DifferentialEquations, Plots, Latexify
+
+struct LaTeXEquation
+    content::String
+end
+
+function Base.show(io::IO, ::MIME"text/latex", x::LaTeXEquation)
+    # Wrap in $$ for display math printing
+    return print(io, "\$\$ " * x.content * " \$\$")
+end
+
+Latexify.set_default(; starred=true)
+
+tspan = (0., 8000.)
+oprob = ODEProblem(kimchi_2020, [], tspan, [])
+sol = solve(oprob)
+```
+
+```@example kimchi_2020
+odesys = convert(ODESystem, kimchi_2020) # hide
+eq = latexify(odesys) # hide
+LaTeXEquation(eq) # hide
+```
+
+Using the parameters specified in the manuscript we obtain the following solution:
+
+```@example kimchi_2020
+plot(sol) # hide
+```

src/BiologicalOscillations.jl

@@ -5,19 +5,21 @@ using Statistics, DSP, Peaks, LatinHypercubeSampling, DataFrames
 using Combinatorics, Luxor, Images, FileIO, Random
 
 # Models
-export elowitz_2000, guan_2008
+export elowitz_2000, guan_2008, kimchi_2020
 # Simulation
 export generate_parameter_sets, equilibrate_ODEs, simulate_ODEs, calculate_simulation_times
 export calculate_oscillatory_status, generate_find_oscillations_output
 # Feature calculation
 export calculate_main_frequency, calculate_amplitude, is_ODE_oscillatory
 # Protein interaction network
 export protein_interaction_network, pin_parameters, pin_timescale, pin_parameter_sets
-export pin_equilibration_times, find_pin_oscillations, pin_nodes_edges
-export pin_hit_rate
+export pin_equilibration_times, find_pin_oscillations, pin_nodes_edges, pin_hit_rate
 # Gene regulatory network
 export gene_regulatory_network, grn_parameters, grn_timescale, grn_parameter_sets
 export grn_equilibration_times, find_grn_oscillations
+# Split enzyme network
+export split_enzyme_network, sen_parameters, sen_timescale, sen_parameter_sets
+export sen_equilibration_times, find_sen_oscillations, sen_hit_rate
 # Network utilities
 export network_permutations, is_same_network, all_network_additions, unique_network_additions
 export is_directed_cycle_graph, is_same_set_of_networks, unique_cycle_addition
@@ -29,6 +31,7 @@ export is_valid_connectivity, connectivity_string_to_matrix
 # Default hyperparameters
 export DEFAULT_PIN_HYPERPARAMETERS, DEFAULT_PIN_PARAMETER_LIMITS, DEFAULT_PIN_SAMPLING_SCALES
 export DEFAULT_GRN_HYPERPARAMETERS, DEFAULT_GRN_PARAMETER_LIMITS, DEFAULT_GRN_SAMPLING_SCALES
+export DEFAULT_SEN_HYPERPARAMETERS, DEFAULT_SEN_PARAMETER_LIMITS, DEFAULT_SEN_SAMPLING_SCALES
 export DEFAULT_SIMULATION_OUTPUT
 # Visualization utilities
 export draw_connectivity
@@ -39,6 +42,7 @@ include("simulation.jl")
 include("network_utilities.jl")
 include("user_input_handling.jl")
 include("feature_calculation.jl")
+include("split_enzyme_network.jl")
 include("gene_regulatory_network.jl")
 include("protein_interaction_network.jl")
 include("draw_utilities.jl")

src/default_hyperparameters.jl

@@ -164,4 +164,59 @@ DEFAULT_GRN_HYPERPARAMETERS = Dict(
     "power_threshold" => 1e-7, # Minimum spectral power that the main peak has to have to be declared oscillatory.
     "amp_variation_threshold" => 0.05, # Maximum tolerated value for peak/trough variation to be declared oscillatory.
     "simulation_output" => DEFAULT_SIMULATION_OUTPUT,
+)
+
+
+"""
+    DEFAULT_SEN_PARAMETER_LIMITS::Dict{String, Tuple{Float64, Float64}}
+
+Default parameter limits for the generation of sen parameter sets. Keys are the names of the parameters, values are tuples of the form (lower_limit, upper_limit).
+"""
+DEFAULT_SEN_PARAMETER_LIMITS = Dict(
+    "k_b" => (1e-2, 1e0), 
+    "k_u" => (1e-3, 1e3), 
+    "n" => (2.0, 4.0), 
+    "κ_tot" => (1e-3, 1e2), 
+    "η" => (1e-1, 1e1), 
+    "P̃" => (1e-3, 1e2)
+)
+
+
+"""
+    DEFAULT_SEN_SAMPLING_SCALES::Dict{String, String}
+
+Default sampling scales for the generation of sen parameter sets. Keys are the names of the parameters, values are strings indicating the sampling scale (log or linear).
+"""
+DEFAULT_SEN_SAMPLING_SCALES = Dict(
+    "k_b" => "log", 
+    "k_u" => "log", 
+    "n" => "linear", 
+    "κ_tot" => "log", 
+    "η" => "log", 
+    "P̃" => "log"
+)
+
+
+"""
+    DEFAULT_SEN_HYPERPARAMETERS
+
+Default hyperparameters for the various algorithms involved in the find_sen_oscillations pipeline.
+"""
+DEFAULT_SEN_HYPERPARAMETERS = Dict(
+    "random_seed" => 1234, # Random seed used to generate parameter sets.
+    "initial_conditions" => NaN, # Set of initial conditions for `find_sen_oscillations`. Size should be equal to the number of samples indicated. If NaN, all species are initialized at 0.5.
+    "parameter_limits" => DEFAULT_SEN_PARAMETER_LIMITS, # Dictionary of parameter limits (lower and upper bound) for each type of parameter in the model.
+    "sampling_scales" => DEFAULT_SEN_SAMPLING_SCALES, # Dictionary of sampling scales for each parameter in the model. Options are "log" or "linear".
+    "sampling_style" => "lhc", # Sampling style. Options are "lhc" (Latin hypercube) or "random".
+    "equilibration_time_multiplier" => 10, # Factor by which to multiply the slowest timescale to get the equilibration time
+    "solver" => RadauIIA5(), # Differential equation solver
+    "abstol" => 1e-7, # Absolute tolerance for the differential equation solver
+    "reltol" => 1e-4, # Relative tolerance for the differential equation solver
+    "maxiters" => 1e7, # Maximum number of iterations for the differential equation solver
+    "simulation_time_multiplier" => 10, # Number of periods to simulate. Used in `calculate_simulation_times`
+    "fft_multiplier" => 100, # Multiplier applied to the number of timepoints in an ODE solution to determine the main frequency. Used in `simulate_ODEs` when calling `calculate_main_frequency`.
+    "freq_variation_threshold" => 0.05, # Maximum tolerated variation in frequency between species in the solution to be declared oscillatory.
+    "power_threshold" => 1e-7, # Minimum spectral power that the main peak has to have to be declared oscillatory.
+    "amp_variation_threshold" => 0.05, # Maximum tolerated value for peak/trough variation to be declared oscillatory.
+    "simulation_output" => DEFAULT_SIMULATION_OUTPUT,
 )

src/models.jl

@@ -39,4 +39,21 @@ guan_2008 = @reaction_network Guan_2008 begin
     (1/sqrt(r))*(a_T + hill(abs(C), b_T, K_T, n_T)), C --> ∅
     # Cdk1 deactivation by Wee1
     sqrt(r)*(a_W + hillr(abs(C), b_W, K_W, n_W)), C --> ∅
+end
+
+
+"""
+    kimchi_2020
+Posttranslational oscillator based on [Kimchi et al. 2020](https://www.science.org/doi/10.1126/sciadv.abc1939)
+"""
+kimchi_2020 = @reaction_network Kimchi_2020 begin
+    @species K(t)=0.5 P(t)=0.5
+    # Parameter values estimated from Fig S2A
+    @parameters p_tilde=1e-9 k_bk=1.0 k_tot=1.0 n=2.0 eta_k=1.0 k_uk=1e-3 k_bp=1e-1 p_tot=7.0 m=2.0 eta_p=1.0 k_up=1e-1
+    # Kinase reactions
+    k_bk * ((k_tot - n * K) / (1 + eta_k * (K / (P + p_tilde)))) ^ n, ∅ --> K
+    k_uk, K --> ∅
+    # Phosphatase reactions
+    k_bp * ((p_tot - m * P) / (1 + eta_p * (K / (P + p_tilde)))) ^ m, ∅ --> P
+    k_up, P --> ∅
 end

src/protein_interaction_network.jl

@@ -297,7 +297,7 @@ end
 
 
 """
-    pin_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target_oscillations::Int=100, max_samples::Int=1000000, max_trials::Int=5, hyperparameters=DEFAULT_PIN_HYPERPARAMETERS)
+    pin_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target_oscillations::Int=100, max_samples::Int=1000000, max_trials::Int=5, hyperparameters=DEFAULT_PIN_HYPERPARAMETERS, verbose=true)
 
 Estimates the hit rate of the oscillatory parameter set search in a protein interaction network
 
@@ -317,9 +317,9 @@ Estimates the hit rate of the oscillatory parameter set search in a protein inte
 """
 function pin_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target_oscillations::Int=100, max_samples::Int=1000000, max_trials::Int=5, hyperparameters=DEFAULT_PIN_HYPERPARAMETERS, verbose=true)
     pin_result = find_pin_oscillations(connectivity, initial_samples, hyperparameters=hyperparameters)
-    oscillatory = sum(pin_result["simulation_result"][!,"is_oscillatory"])
+    oscillatory = sum(pin_result["simulation_result"][!, "is_oscillatory"])
     if verbose
-        println("Oscillatory: $oscillatory, Samples $initial_samples")
+        println("Oscillatory: $oscillatory, Samples: $initial_samples")
     end
     hit_rate =  oscillatory / initial_samples
 
@@ -337,7 +337,7 @@ function pin_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target
             pin_result = find_pin_oscillations(connectivity, samples, hyperparameters=hyperparameters)
             oscillatory = sum(pin_result["simulation_result"][!,"is_oscillatory"])
             if verbose
-                println("Oscillatory: $oscillatory, Samples $samples")
+                println("Oscillatory: $oscillatory, Samples: $samples")
             end
             hit_rate = oscillatory / samples
             trial += 1

src/simulation.jl

@@ -29,7 +29,7 @@ function generate_parameter_sets(samples::Int, parameter_limits::AbstractVector,
     scaling_plan = Tuple{Float64,Float64}[]
     for (idx, limits) in enumerate(parameter_limits)
         if sampling_scales[idx] == "linear"
-            scaling_plan = vcat(scaling_plan, [limits])
+            scaling_plan = vcat(scaling_plan, [Tuple(limits)])
         elseif sampling_scales[idx] == "log"
             scaling_plan = vcat(scaling_plan, [(log10(limits[1]), log10(limits[2]))])
         end