From 50f036f647f2a05a45f2be88583a5cc0617538fc Mon Sep 17 00:00:00 2001 From: ftavella Date: Tue, 9 Jan 2024 11:13:04 -0500 Subject: [PATCH 1/6] First implementation of two-equation Kimchi et al. 2020 model --- docs/pages.jl | 2 +- docs/src/models/kimchi_2020.md | 34 ++++++++++++++++++++++++++++++++++ src/BiologicalOscillations.jl | 2 +- src/models.jl | 17 +++++++++++++++++ 4 files changed, 53 insertions(+), 2 deletions(-) create mode 100644 docs/src/models/kimchi_2020.md diff --git a/docs/pages.jl b/docs/pages.jl index fe0b0e2..703452b 100644 --- a/docs/pages.jl +++ b/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"], ] \ No newline at end of file diff --git a/docs/src/models/kimchi_2020.md b/docs/src/models/kimchi_2020.md new file mode 100644 index 0000000..4788774 --- /dev/null +++ b/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 +``` \ No newline at end of file diff --git a/src/BiologicalOscillations.jl b/src/BiologicalOscillations.jl index 39fb39d..277197e 100644 --- a/src/BiologicalOscillations.jl +++ b/src/BiologicalOscillations.jl @@ -5,7 +5,7 @@ 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 diff --git a/src/models.jl b/src/models.jl index 37aed30..5984e6a 100644 --- a/src/models.jl +++ b/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 \ No newline at end of file From c56494d53393165d135941cf2323de6ae93aad74 Mon Sep 17 00:00:00 2001 From: ftavella Date: Thu, 11 Jan 2024 16:44:26 -0500 Subject: [PATCH 2/6] Split enzyme network implementation plus minor bug fixes --- src/BiologicalOscillations.jl | 8 +- src/default_hyperparameters.jl | 55 +++ src/protein_interaction_network.jl | 8 +- src/simulation.jl | 2 +- src/split_enzyme_network.jl | 363 ++++++++++++++++++ test/protein_interaction_network/pin_tests.jl | 6 +- test/runtests.jl | 1 + test/split_enzyme_network/sen_test.jl | 358 +++++++++++++++++ 8 files changed, 791 insertions(+), 10 deletions(-) create mode 100644 src/split_enzyme_network.jl create mode 100644 test/split_enzyme_network/sen_test.jl diff --git a/src/BiologicalOscillations.jl b/src/BiologicalOscillations.jl index 277197e..24051bf 100644 --- a/src/BiologicalOscillations.jl +++ b/src/BiologicalOscillations.jl @@ -13,11 +13,13 @@ export calculate_oscillatory_status, generate_find_oscillations_output 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") diff --git a/src/default_hyperparameters.jl b/src/default_hyperparameters.jl index 5065a1f..1d33d46 100644 --- a/src/default_hyperparameters.jl +++ b/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, 1e-3), + "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, ) \ No newline at end of file diff --git a/src/protein_interaction_network.jl b/src/protein_interaction_network.jl index 5dca1d5..8f6c34f 100644 --- a/src/protein_interaction_network.jl +++ b/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 diff --git a/src/simulation.jl b/src/simulation.jl index 9636e72..aa3d30f 100644 --- a/src/simulation.jl +++ b/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 diff --git a/src/split_enzyme_network.jl b/src/split_enzyme_network.jl new file mode 100644 index 0000000..982aa20 --- /dev/null +++ b/src/split_enzyme_network.jl @@ -0,0 +1,363 @@ +function split_enzyme_monomer_concentration(E, K, P, monomer_tot, n_monomers, η) + """ + Concentration of enzyme monomers under the assumption of Michaelis-Menten dynamics plus a separation of timescales with + phosphorylation/dephosphorylation reactions equilibrating much faster than split enzyme self-assembly. + Expression taken from [Kimchi et al. 2020](https://www.science.org/doi/10.1126/sciadv.abc1939) + + # Arguments + - `E`: Concentration of full enzyme + - `K`: Concentration of kinase sequestering monomers + - `P`: Concentration of phosphatase freeing monomers + - `monomer_tot`: Total concentration of enzyme monomers + - `n_monomers`: Number of monomers in a single split enzyme + - `η`: Ratio of specificities between kinase and phosphatase + """ + return abs((monomer_tot - n_monomers * E) / (1 + η * (K / P))) +end + + +""" + nodes_have_valid_split_enzyme_outputs(connectivity::AbstractMatrix) + +Check if each node has only one type of output (kinase or phosphatase). This translates to checking that each column of the connectivity matrix has outputs of the same sign. + +# Arguments +- `connectivity::AbstractMatrix`: A 2-dimensional matrix filled with -1, 0, and 1 values indicating the edges of the network. + -1 indicates a kinase input, 1 indicates a phosphatase input, and 0 indicates no input. Please note that a node can only be + a kinase or a phosphatase, not both. +""" +function nodes_have_valid_split_enzyme_outputs(connectivity::AbstractMatrix) + result = true + errmsg = "" + # Check that each column has the same sign + for i in axes(connectivity, 2) + column = connectivity[:,i] + if !(all(column .>= 0) || all(column .<= 0)) + result = false + errmsg = "Invalid node output found. Each column of the connectivity matrix must have the same sign" + return [result, errmsg] + end + end + return [result, errmsg] +end + + +""" + needs_constitutive_kinase(connectivity::AbstractMatrix) + +Check if there is at least one node that doesn't have a kinase input (i.e. a node that doesn't have a -1 in its row). + +# Arguments +- `connectivity::AbstractMatrix`: A 2 dimensional matrix filled with -1, 0, and 1 values indicating the edges of the network. + -1 indicates a kinase input, 1 indicates a phosphatase input, and 0 indicates no input. +""" +function needs_constitutive_kinase(connectivity::AbstractMatrix) + needs = false + for i in axes(connectivity, 1) + row = connectivity[i,:] + if !any(row .== -1) + needs = true + end + end + return needs +end + + +""" + split_enzyme_network(connectivity::AbstractMatrix) + +Creates a ReactionSystem of interacting split enzymes based on the connectivity matrix. +The model is based on [Kimchi et al. 2020](https://www.science.org/doi/10.1126/sciadv.abc1939). +A constitutive phosphatase parameter (P̃) is added to all nodes. A constitutive kinase (K̃) is +only added to those nodes that don't have a kinase input. + +# Arguments +- `connectivity::AbstractMatrix`: A 2 dimensional matrix filled with -1, 0, and 1 values indicating the edges of the network. + -1 indicates a kinase input, 1 indicates a phosphatase input, and 0 indicates no input. Please note that a node can only be + a kinase or a phosphatase, not both. This means that all the non-zero values of a column must have the same sign. +""" +function split_enzyme_network(connectivity::AbstractMatrix) + # Check that the input is correct + is_valid, errmsg = is_valid_connectivity(connectivity) + has_correct_node_outputs, node_output_errmsg = nodes_have_valid_split_enzyme_outputs(connectivity) + if !is_valid + throw(DomainError(connectivity, errmsg)) + elseif !has_correct_node_outputs + throw(DomainError(connectivity, node_output_errmsg)) + end + # Number of nodes + N = size(connectivity, 1) + # Does the network need a constitutive kinase? + constitutive_kinase_needed = needs_constitutive_kinase(connectivity) + @variables t + @species (X(t))[collect(1:N)] + if constitutive_kinase_needed + @parameters k_b[1:N], k_u[1:N], n[1:N], κ_tot[1:N], η, K̃, P̃ + else + @parameters k_b[1:N], k_u[1:N], n[1:N], κ_tot[1:N], η, P̃ + end + rxs = Reaction[] + for i=1:N + # Search for all inputs to node i + K = 0 + P = P̃ + for j=1:N + if connectivity[i,j] == -1 + K += X[j] + elseif connectivity[i,j] == 1 + P += X[j] + end + end + if iszero(K) && constitutive_kinase_needed + K += K̃ + end + # Define binding and unbinding reactions + monomer = split_enzyme_monomer_concentration(X[i], K, P, κ_tot[i], n[i], η) + push!(rxs, Reaction(k_b[i] * monomer ^ n[i], nothing, [X[i]])) + push!(rxs, Reaction(k_u[i], [X[i]], nothing)) + end + # Reorder parameters + parameters = [[k_b[i], k_u[i], n[i], κ_tot[i]] for i=1:N] + parameters = vcat(parameters...) + parameters = vcat(parameters, η, P̃) + if constitutive_kinase_needed + parameters = vcat(parameters, K̃) + end + @named model = ReactionSystem(rxs, t, X, parameters) + return model +end + + +""" + sen_parameter_sets(model::ReactionSystem, samples::Int, random_seed::Int; parameter_limits=DEFAULT_SEN_PARAMETER_LIMITS, sampling_scales=DEFAULT_SEN_SAMPLING_SCALES, sampling_style="lhc") + +Creates an array of paramter sets for a split enzyme network model. + +# Arguments (Required) +- `model::ReactionSystem`: A split enzyme network model +- `samples::Int`: Number of parameter sets to generate +- `random_seed::Int`: Random seed for the sampling algorithm reproducibility + +# Arguments (Optional) +- `parameter_limits::Dict`: Dictionary of parameter limits (lower and upper bound) for each type of parameter in the model. +- `sampling_scales::Dict`: Dictionary of sampling scales for each parameter in the model. Options are "log" or "linear". +- `sampling_style::String`: Sampling style. Options are "lhc" (Latin hypercube) or "random". + +""" +function sen_parameter_sets(model::ReactionSystem, samples::Int, random_seed::Int; parameter_limits=DEFAULT_SEN_PARAMETER_LIMITS, sampling_scales=DEFAULT_SEN_SAMPLING_SCALES, sampling_style="lhc") + N = length(species(model)) + P = length(parameters(model)) + + limits = [] + for i in 1:N + push!(limits, parameter_limits["k_b"]) + push!(limits, parameter_limits["k_u"]) + push!(limits, parameter_limits["n"]) + push!(limits, parameter_limits["κ_tot"]) + end + + push!(limits, parameter_limits["η"]) + + for i in 1:(P-4*N-1) + push!(limits, parameter_limits["P̃"]) + end + + scales = [] + for i in 1:N + push!(scales, sampling_scales["k_b"]) + push!(scales, sampling_scales["k_u"]) + push!(scales, sampling_scales["n"]) + push!(scales, sampling_scales["κ_tot"]) + end + + push!(scales, sampling_scales["η"]) + + for i in 1:(P-4*N-1) + push!(scales, sampling_scales["P̃"]) + end + + parameter_array = generate_parameter_sets(samples, limits, scales, random_seed;sampling_style=sampling_style) + + # Force number of monomers to be integer + parameter_array[:,3:4:4*N] = round.(parameter_array[:,3:4:4*N]) + + return parameter_array +end + + +""" + sen_timescale(k_u::AbstractVector, k_b::AbstractVector, n::AbstractVector, κ_tot::AbstractVector) + +Calculate the slowest timescale of a split enzyme network model. + +# Arguments +- `k_u::AbstractVector`: Vector of unbinding rate constants +- `k_b::AbstractVector`: Vector of binding rate constants +- `n::AbstractVector`: Vector of number of monomers in each split enzyme +- `κ_tot::AbstractVector`: Vector of total concentration of enzyme monomers + +# Returns +- `timescale::Real` +""" +function sen_timescale(k_u::AbstractVector, k_b::AbstractVector, n::AbstractVector, κ_tot::AbstractVector) + rates = vcat(k_u, k_b .* (κ_tot) .^ (n .- 1)) + return 1.0 / minimum(rates) +end + + +""" + sen_equilibration_times(model::ReactionSystem, parameter_sets::AbstractArray; equilibration_time_multiplier=10) + +Calculates the equilibration time for each dataset + +# Arguments (Required) +- `model::ReactionSystem`: A split enzyme network model +- `parameter_sets::AbstractArray`: Array of parameter sets generated with [`sen_parameter_sets`](@ref) + +# Arguments (Optional) +- `equilibration_time_multiplier::Real`: Factor by which to multiply the slowest timescale to get the equilibration time + +# Returns +- `equilibration_times::Array{Float64}`: Vector of equilibration times for each parameter set +""" +function sen_equilibration_times(model::ReactionSystem, parameter_sets::AbstractArray; equilibration_time_multiplier=10) + N = length(species(model)) + equilibration_times = [] + for i in axes(parameter_sets, 1) + k_b = parameter_sets[i, 1:4:4*N] + k_u = parameter_sets[i, 2:4:4*N] + n = parameter_sets[i, 3:4:4*N] + κ_tot = parameter_sets[i, 4:4:4*N] + timescale = sen_timescale(k_u, k_b, n, κ_tot) + push!(equilibration_times, timescale * equilibration_time_multiplier) + end + return equilibration_times +end + + +""" + find_sen_oscillations(model::ReactionSystem, samples::Int; hyperparameters=DEFAULT_SEN_HYPERPARAMETERS) +""" +function find_sen_oscillations(connectivity::AbstractMatrix, samples::Int; hyperparameters=DEFAULT_SEN_HYPERPARAMETERS) + # Unpack hyperparameters + random_seed = hyperparameters["random_seed"] + initial_conditions = hyperparameters["initial_conditions"] + parameter_limits = hyperparameters["parameter_limits"] + sampling_scales = hyperparameters["sampling_scales"] + sampling_style = hyperparameters["sampling_style"] + equilibration_time_multiplier = hyperparameters["equilibration_time_multiplier"] + solver = hyperparameters["solver"] + abstol = hyperparameters["abstol"] + reltol = hyperparameters["reltol"] + maxiters = hyperparameters["maxiters"] + simulation_time_multiplier = hyperparameters["simulation_time_multiplier"] + fft_multiplier = hyperparameters["fft_multiplier"] + freq_variation_threshold = hyperparameters["freq_variation_threshold"] + power_threshold = hyperparameters["power_threshold"] + amp_variation_threshold = hyperparameters["amp_variation_threshold"] + + model = split_enzyme_network(connectivity) + N = length(species(model)) + parameter_sets = sen_parameter_sets(model, samples, random_seed; + parameter_limits=parameter_limits, + sampling_scales=sampling_scales, + sampling_style=sampling_style) + equilibration_times = sen_equilibration_times(model, parameter_sets; + equilibration_time_multiplier=equilibration_time_multiplier) + + if isnan(initial_conditions) + initial_conditions = [zeros(N) for i in 1:samples] + end + + # Equilibrate + equilibration_data = equilibrate_ODEs(model, parameter_sets, initial_conditions, equilibration_times; + solver=solver, abstol=abstol, reltol=reltol, maxiters=maxiters) + # Filter out solutions with velocity vector smaller than the mean velocity obtained from the equilibration + velocity = equilibration_data["final_velocity"] + cutoff = 10 ^ mean(log10.(velocity)) + filter = velocity .> cutoff + simulation_times = calculate_simulation_times(equilibration_data, equilibration_times; + simulation_time_multiplier=simulation_time_multiplier) + # Simulate + simulation_data = simulate_ODEs(model, parameter_sets[filter,:], equilibration_data["final_state"][filter], simulation_times[filter]; + solver=solver, abstol=abstol, reltol=reltol, maxiters=maxiters, fft_multiplier=fft_multiplier) + # Check for oscillations + oscillatory_status = calculate_oscillatory_status(simulation_data; + freq_variation_threshold=freq_variation_threshold, + power_threshold=power_threshold, + amp_variation_threshold=amp_variation_threshold) + # Output + result = generate_find_oscillations_output(model, parameter_sets, equilibration_data, equilibration_times, + simulation_data, simulation_times, oscillatory_status, hyperparameters) + return result +end + + +""" + sen_hit_rate(connectivity::AbstractMAtrix, initial_samples::Int; target_oscillations::Int=100, max_samples::Int=1000000, max_trials::Int=5, hyperparameters=DEFAULT_SEN_HYPERPARAMETERS, verbose=true) + +Estimates the hit rate of the oscillatory parameter set search in a split enzyme network model. + +# Arguments (Required) +- `connectivity::AbstractMatrix`: A 2 dimensional matrix filled with -1, 0, and 1 values indicating the edges of the network. + -1 indicates a kinase input, 1 indicates a phosphatase input, and 0 indicates no input. Please note that a node can only be + a kinase or a phosphatase, not both. This means that all the non-zero values of a column must have the same sign. +- `initial_samples::Int`: Number of initial parameter sets to generate + +# Arguments (Optional) +- `target_oscillations::Int`: Number of oscillatory parameter sets to find +- `max_samples::Int`: Maximum number of samples allowed +- `max_trials::Int`: Maximum number of trials allowed +- `hyperparameters::Dict`: Dictionary of hyperparameters for the oscillatory parameter set search +- `verbose::Bool`: Whether to print the progress of the search +""" +function sen_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target_oscillations::Int=100, max_samples::Int=1000000, max_trials::Int=5, hyperparameters=DEFAULT_SEN_HYPERPARAMETERS, verbose=true) + sen_result = find_sen_oscillations(connectivity, initial_samples; hyperparameters=hyperparameters) + oscillatory = sum(sen_result["simulation_result"][!, "is_oscillatory"]) + if verbose + println("Oscillatory: $oscillatory, Samples: $initial_samples") + end + hit_rate = oscillatory / initial_samples + + if oscillatory > target_oscillations + return hit_rate + else + samples = initial_samples + trial = 1 + while oscillatory < target_oscillations && trial < max_trials && samples < max_samples + if oscillatory == 0 + samples *= 2 + else + samples = ceil(Int, samples * (target_oscillations + 10) / oscillatory) + end + sen_result = find_sen_oscillations(connectivity, samples; hyperparameters=hyperparameters) + oscillatory = sum(sen_result["simulation_result"][!, "is_oscillatory"]) + if verbose + println("Oscillatory: $oscillatory, Samples: $samples") + end + hit_rate = oscillatory / samples + trial += 1 + end + return hit_rate + end +end + + +""" +Goodwin-type negative feedback oscillator comprised of two split phosphatases, one split kinase, +one constitutive kinase, and one constitutive phosphatase. Model based on [Kimchi et al. 2020](https://www.science.org/doi/10.1126/sciadv.abc1939) +Connectivity: K phosphorylates P₁, P₁ dephosphorylates P₂, P₂ dephosphorylates K +""" +self_assembly_negative_feedback = @reaction_network Self_assembly_negative_feedback begin + @species K(t)=0.5 P₁(t)=0.5 P₂(t)=0.5 + @parameters k_bκ=1e-1 k_uκ=1e0 n=2.0 κ_tot=3.0 k_bρ₁=1e0 k_uρ₁=1e-3 m=2.0 ρ₁_tot=7.0 k_bρ₂=1e-1 k_uρ₂=1e-1 l=2.0 ρ₂_tot=1.0 η=1.0 P̃=1e-5 K̃=1e-6 + # K reactions + k_bκ * split_enzyme_monomer_concentration(K, K̃, P₂ + P̃, κ_tot, n, η) ^ n, ∅ --> K + k_uκ, K --> ∅ + # P₁ reactions + k_bρ₁ * split_enzyme_monomer_concentration(P₁, K, P̃, ρ₁_tot, m, η) ^ m, ∅ --> P₁ + k_uρ₁, P₁ --> ∅ + # P₂ reactions + k_bρ₂ * split_enzyme_monomer_concentration(P₂, K̃, P₁ + P̃, ρ₂_tot, l, η) ^ l, ∅ --> P₂ + k_uρ₂, P₂ --> ∅ +end diff --git a/test/protein_interaction_network/pin_tests.jl b/test/protein_interaction_network/pin_tests.jl index 1d1bc59..652e253 100644 --- a/test/protein_interaction_network/pin_tests.jl +++ b/test/protein_interaction_network/pin_tests.jl @@ -162,7 +162,7 @@ calculated_hit_rate = pin_hit_rate(connectivity_P0_6, samples; verbose=false) samples = 200 random_seed = 4567 connectivity_T0 = [0 0 -1;-1 0 0;0 -1 0] -hyperparameters = DEFAULT_PIN_HYPERPARAMETERS +hyperparameters = deepcopy(DEFAULT_PIN_HYPERPARAMETERS) hyperparameters["random_seed"] = random_seed result = find_pin_oscillations(connectivity_T0, samples; hyperparameters=hyperparameters) @@ -200,7 +200,7 @@ end samples = 200 random_seed = 4567 connectivity_T0 = [0 0 -1;-1 0 0;0 -1 0] -hyperparameters = DEFAULT_PIN_HYPERPARAMETERS +hyperparameters = deepcopy(DEFAULT_PIN_HYPERPARAMETERS) hyperparameters["random_seed"] = random_seed solver = hyperparameters["solver"] result = find_pin_oscillations(connectivity_T0, samples; @@ -250,7 +250,7 @@ end # Test that output can be customized samples = 200 connectivity_T0 = [0 0 -1;-1 0 0;0 -1 0] -hyperparameters = DEFAULT_PIN_HYPERPARAMETERS +hyperparameters = deepcopy(DEFAULT_PIN_HYPERPARAMETERS) sim_output_config = Dict( "model" => true, diff --git a/test/runtests.jl b/test/runtests.jl index e1377db..1afa86b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,7 @@ using Test @testset "Protein Interaction Networks" begin include("protein_interaction_network/pin_tests.jl") end @testset "Gene Regulatory Networks" begin include("gene_regulatory_network/grn_tests.jl") end +@testset "Split Enzyme Networks" begin include("split_enzyme_network/sen_tests.jl") end @testset "User Input Handling" begin include("user_input_handling_tests.jl") end @testset "Feature Calculation" begin include("feature_calculation_tests.jl") end @testset "Network Utilities" begin include("network_utilities_tests.jl") end diff --git a/test/split_enzyme_network/sen_test.jl b/test/split_enzyme_network/sen_test.jl new file mode 100644 index 0000000..87f253f --- /dev/null +++ b/test/split_enzyme_network/sen_test.jl @@ -0,0 +1,358 @@ +using BiologicalOscillations, Catalyst, ModelingToolkit, DifferentialEquations +using DataFrames, Peaks + +# Test `nodes_have_valid_split_enzyme_outputs` +invalid_connectivity_1 = [1 0 1; 0 1 0; 0 0 -1] +invalid_connectivity_2 = [1 -1 -1; 0 1 0; 0 0 -1] +invalid_connectivity_3 = [1 1 -1 0; 0 1 0 1; 0 0 -1 1;0 0 0 -1] +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(invalid_connectivity_1) +@test result == false +@test errmsg == "Invalid node output found. Each column of the connectivity matrix must have the same sign" +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(invalid_connectivity_2) +@test result == false +@test errmsg == "Invalid node output found. Each column of the connectivity matrix must have the same sign" +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(invalid_connectivity_3) +@test result == false +@test errmsg == "Invalid node output found. Each column of the connectivity matrix must have the same sign" + +valid_connectivity_1 = [1 0 1; 0 1 0; 0 0 1] +valid_connectivity_2 = [1 1 -1; 0 1 0; 0 0 -1] +valid_connectivity_3 = [1 1 -1 0; 0 1 0 1; 0 0 -1 1;0 0 0 1] +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(valid_connectivity_1) +@test result == true +@test errmsg == "" +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(valid_connectivity_2) +@test result == true +@test errmsg == "" +result, errmsg = BiologicalOscillations.nodes_have_valid_split_enzyme_outputs(valid_connectivity_3) +@test result == true +@test errmsg == "" + +# Test `needs_constitutive_kinase` +needs_kinase_1 = [0 0 1; 1 0 0; 0 1 0] +needs_kinase_2 = [0 0 1; -1 0 0; 0 -1 0] +needs_kinase_3 = [0 0 1; 1 0 0; 0 -1 0] +result = BiologicalOscillations.needs_constitutive_kinase(needs_kinase_1) +@test result == true +result = BiologicalOscillations.needs_constitutive_kinase(needs_kinase_2) +@test result == true +result = BiologicalOscillations.needs_constitutive_kinase(needs_kinase_3) +@test result == true + +doesnt_need_kinase_1 = [0 0 -1; -1 0 0; 0 -1 0] +doesnt_need_kinase_2 = [1 0 -1; 1 -1 0; 0 -1 0] +doesnt_need_kinase_3 = [1 0 -1; 1 -1 0; 1 -1 -1] +result = BiologicalOscillations.needs_constitutive_kinase(doesnt_need_kinase_1) +@test result == false +result = BiologicalOscillations.needs_constitutive_kinase(doesnt_need_kinase_2) +@test result == false +result = BiologicalOscillations.needs_constitutive_kinase(doesnt_need_kinase_3) +@test result == false + +# Test `split_enzyme_network` +## Error handling +@test_throws MethodError split_enzyme_network([]) +@test_throws MethodError split_enzyme_network("ASD") +@test_throws DomainError split_enzyme_network(zeros((0,0))) +@test_throws DomainError split_enzyme_network([0 1 -1]) +@test_throws DomainError split_enzyme_network(ones((1,2))) +@test_throws DomainError split_enzyme_network([0 1;3 -1]) +## Correct reactions - Repressilator +@variables t +@species (X(t))[collect(1:3)] +@parameters k_b[1:3], k_u[1:3], n[1:3], κ_tot[1:3], η, P̃ + +true_reactions = [] +for i=1:3 + K = X[mod(i-2,3)+1] # Get the previous species in the loop + monomer = BiologicalOscillations.split_enzyme_monomer_concentration(X[i], K, P̃, κ_tot[i], n[i], η) + push!(true_reactions, Reaction(k_b[i] * monomer ^ n[i], nothing, [X[i]])) + push!(true_reactions, Reaction(k_u[i], [X[i]], nothing)) +end +@named true_repressilator = ReactionSystem(true_reactions, t) +generated_repressilator = split_enzyme_network([0 0 -1;-1 0 0;0 -1 0]) +@test reactions(generated_repressilator) == reactions(true_repressilator) +## Correct reactions - Goodwin + Incoherent positive feedback loop +@parameters k_b[1:3], k_u[1:3], n[1:3], κ_tot[1:3], η, P̃, K̃ +true_reactions = [] +monomer_1 = BiologicalOscillations.split_enzyme_monomer_concentration(X[1], K̃, X[3] + P̃, κ_tot[1], n[1], η) +push!(true_reactions, Reaction(k_b[1] * monomer_1 ^ n[1], nothing, [X[1]])) +push!(true_reactions, Reaction(k_u[1], [X[1]], nothing)) +monomer_2 = BiologicalOscillations.split_enzyme_monomer_concentration(X[2], K̃, X[1] + P̃, κ_tot[2], n[2], η) +push!(true_reactions, Reaction(k_b[2] * monomer_2 ^ n[2], nothing, [X[2]])) +push!(true_reactions, Reaction(k_u[2], [X[2]], nothing)) +monomer_3 = BiologicalOscillations.split_enzyme_monomer_concentration(X[3], X[2], X[3] + P̃, κ_tot[3], n[3], η) +push!(true_reactions, Reaction(k_b[3] * monomer_3 ^ n[3], nothing, [X[3]])) +push!(true_reactions, Reaction(k_u[3], [X[3]], nothing)) + +@named true_goodwin = ReactionSystem(true_reactions, t) +generated_goodwin = split_enzyme_network([0 0 1;1 0 0;0 -1 1]) +@test reactions(generated_goodwin) == reactions(true_goodwin) + +## Test `sen_timescale` +k_u = [1, 2, 3] +k_b = [4, 5, 6] +n = [7, 8, 9] +κ_tot = [10, 11, 12] +timescale = BiologicalOscillations.sen_timescale(k_u, k_b, n, κ_tot) +@test timescale == 1.0 + +## Test `sen_parameter_sets` +repressilator = split_enzyme_network([0 0 -1;-1 0 0;0 -1 0]) +samples = 10 +parameter_array = sen_parameter_sets(repressilator, samples, 123) +n_parameters = length(parameters(repressilator)) +@test size(parameter_array) == (samples, n_parameters) +n_values = parameter_array[:, 3:4:4*3] +@test all(n_values .== round.(n_values)) # Check that n is an integer + +## Test `sen_equilibration_times` +repressilator = split_enzyme_network([0 0 -1;-1 0 0;0 -1 0]) +samples = 10 +parameter_array = sen_parameter_sets(repressilator, samples, 123) +equilibration_times = sen_equilibration_times(repressilator, parameter_array) +@test size(equilibration_times) == (samples,) +# TODO: Create a known solution and test that the equilibration times are correct + +## Test `find_sen_oscillations` +# TODO: Postponing until runtime is reduced +""" +samples = 5000 +connectivity_T2 = [0 0 1;1 0 0;0 -1 0] +T2_hit_rate = 0.0012 +sen_result_T2 = find_sen_oscillations(connectivity_T2, samples) +oscillatory_df = filter(row -> row["is_oscillatory"] == true, sen_result_T2["simulation_result"]) +oscillatory_solutions = size(oscillatory_df, 1) +@test oscillatory_solutions > samples*T2_hit_rate*0.8 + +samples = 5000 +connectivity_T2_12 = [0 0 1; 1 0 0; 0 -1 1] +T2_12_hit_rate = 0.0024 +sen_result_T2_12 = find_sen_oscillations(connectivity_T2_12, samples) +oscillatory_df = filter(row -> row["is_oscillatory"] == true, sen_result_T2_12["simulation_result"]) +oscillatory_solutions = size(oscillatory_df, 1) +@test oscillatory_solutions > samples*T2_12_hit_rate*0.8 + +samples = 1000 +connectivity_P0 = [0 0 0 0 -1; -1 0 0 0 0; 0 -1 0 0 0; 0 0 -1 0 0; 0 0 0 -1 0] +P0_hit_rate = 0.04 +sen_result_P0 = find_sen_oscillations(connectivity_P0, samples) +oscillatory_df = filter(row -> row["is_oscillatory"] == true, sen_result_P0["simulation_result"]) +oscillatory_solutions = size(oscillatory_df, 1) +@test oscillatory_solutions > samples*P0_hit_rate*0.8 +""" + +## Test `sen_hit_rate` +# TODO: Postponing until runtime is reduced and simulations validated +""" +samples = 5000 +connectivity_T2 = [0 0 1;1 0 0;0 -1 0] +T2_hit_rate = 0.0012 +calculated_hit_rate = sen_hit_rate(connectivity_T2, samples; verbose=false) +@test calculated_hit_rate ≈ T2_hit_rate rtol=0.3 + +samples = 5000 +connectivity_T2_12 = [0 0 1; 1 0 0; 0 -1 1] +T2_12_hit_rate = 0.0024 +calculated_hit_rate = sen_hit_rate(connectivity_T2_12, samples; verbose=false) +@test calculated_hit_rate ≈ T2_12_hit_rate rtol=0.3 + +samples = 1000 +connectivity_P0 = [0 0 0 0 -1; -1 0 0 0 0; 0 -1 0 0 0; 0 0 -1 0 0; 0 0 0 -1 0] +P0_hit_rate = 0.04 +calculated_hit_rate = sen_hit_rate(connectivity_P0, samples; verbose=false) +@test calculated_hit_rate ≈ P0_hit_rate rtol=0.3 +""" + +## Test that parameters from run and generated with a random seed are the same +samples = 200 +random_seed = 5678 +connectivity_T2 = [0 0 1;1 0 0;0 -1 0] +hyperparameters = deepcopy(DEFAULT_SEN_HYPERPARAMETERS) +hyperparameters["random_seed"] = random_seed +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) + +oscillatory_params = result["parameter_sets"]["oscillatory"] +fixed_point_params = result["parameter_sets"]["non_oscillatory"] +all_params = vcat(oscillatory_params, fixed_point_params) +all_params = sort(all_params, :parameter_index) +all_params = all_params[:, 2:end] +column_order = ["k_b[1]", "k_u[1]", "n[1]", "κ_tot[1]", + "k_b[2]", "k_u[2]", "n[2]", "κ_tot[2]", + "k_b[3]", "k_u[3]", "n[3]", "κ_tot[3]", + "η", "P̃", "K̃"] +all_params = select(all_params, column_order) + +model = result["model"] +parameter_limits = hyperparameters["parameter_limits"] +sampling_scales = hyperparameters["sampling_scales"] +sampling_style = hyperparameters["sampling_style"] + +parameter_sets = sen_parameter_sets(model, samples, random_seed; + parameter_limits=parameter_limits, + sampling_scales=sampling_scales, + sampling_style=sampling_style) + +for (idx, row) in enumerate(eachrow(all_params)) + params_1 = parameter_sets[idx, :] + params_2 = collect(row) + @test params_1 == params_2 +end + +# Test that parameters are truly oscillatory and non-oscillatory +# Values taken from run on 01/11/2024 +samples = 1500 +random_seed = 9283 +connectivity_T2_12 = [0 0 1; 1 0 0; 0 -1 1] +hyperparameters = deepcopy(DEFAULT_SEN_HYPERPARAMETERS) +hyperparameters["random_seed"] = random_seed +solver = hyperparameters["solver"] +result = find_sen_oscillations(connectivity_T2_12, samples; + hyperparameters=hyperparameters) +model = result["model"] +oscillatory_params = result["parameter_sets"]["oscillatory"] + +peak_number_found = zeros(5) + +for selected_idx in 1:5 + parameter_index = oscillatory_params[selected_idx, 1] + final_states = select(result["simulation_result"], r"parameter_index|final_state_") + final_states = filter(row -> row["parameter_index"] == parameter_index, final_states) + initial_condition = collect(final_states[1, 1:end-1]) + final_time = 50000 + params = collect(oscillatory_params[selected_idx, 2:end]) + ode_problem = ODEProblem(model, initial_condition, + (0.0, final_time), params) + sol = solve(ode_problem, solver, saveat = final_time/5000) + df = DataFrame(t = sol.t, x = sol[1, :], y = sol[2, :], z = sol[3, :]) + peak_num = size(findmaxima(df.x)[1], 1) + peak_number_found[selected_idx] = peak_num +end + +expected_peaks = [51, 5, 569, 14, 3] +@test peak_number_found == expected_peaks +# Test that parameters are truly non-oscillatory +non_oscillatory_params = result["parameter_sets"]["non_oscillatory"] + +for selected_idx in 1:5 + parameter_index = non_oscillatory_params[selected_idx, 1] + final_states = select(result["equilibration_result"], r"parameter_index|final_state_") + final_states = filter(row -> row["parameter_index"] == parameter_index, final_states) + final_condition = collect(final_states[1, 1:end-1]) + initial_condition = zeros(3) + final_time = 50000 + parameters = collect(non_oscillatory_params[selected_idx, 2:end]) + ode_problem = ODEProblem(model, initial_condition, + (0.0, final_time), parameters) + sol = solve(ode_problem, solver, saveat = final_time/5000) + df = DataFrame(t = sol.t, x = sol[1, :], y = sol[2, :], z = sol[3, :]) + final_state_simulation = collect(df[end, 2:end]) + @test final_state_simulation ≈ final_condition rtol=1e-4 +end + +# Test that the output can be customized +samples = 200 +connectivity_T2 = [0 0 1;1 0 0;0 -1 0] +hyperparameters = deepcopy(DEFAULT_SEN_HYPERPARAMETERS) + +sim_output_config = Dict( + "model" => true, +) +hyperparameters["simulation_output"] = sim_output_config +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) +@test result["model"] isa ReactionSystem + +sim_output_config = Dict( + "hyperparameters" => true, +) +hyperparameters["simulation_output"] = sim_output_config +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) +@test result["hyperparameters"] isa Dict +@test keys(result["hyperparameters"]) == keys(hyperparameters) +for key in keys(hyperparameters) + # don't compare initial_conditions because it's NaN + if key == "initial_conditions" + continue + end + @test result["hyperparameters"][key] == hyperparameters[key] +end + +sim_output_config = Dict( + "parameter_sets" => Dict( + "oscillatory" => true, + "non_oscillatory" => true + ) +) +hyperparameters["simulation_output"] = sim_output_config +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) +@test result["parameter_sets"] isa Dict +@test result["parameter_sets"]["oscillatory"] isa DataFrame +@test result["parameter_sets"]["non_oscillatory"] isa DataFrame + +sim_output_config = Dict( + "equilibration_result" => Dict( + "parameter_index" => true, + "equilibration_time" => true, + "final_velocity" => true, + "final_state" => true, + "frequency" => true, + "is_steady_state" => true + ) +) +hyperparameters["simulation_output"] = sim_output_config +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) +@test result["equilibration_result"] isa DataFrame +@test size(result["equilibration_result"], 2) == 8 + +sim_output_config = Dict( + "simulation_result" => Dict( + "parameter_index" => true, + "simulation_time" => true, + "final_state" => true, + "frequency" => true, + "fft_power" => true, + "amplitude" => true, + "peak_variation" => true, + "trough_variation" => true, + "is_oscillatory" => true + ) +) + +hyperparameters["simulation_output"] = sim_output_config +result = find_sen_oscillations(connectivity_T2, samples; + hyperparameters=hyperparameters) +@test result["simulation_result"] isa DataFrame +@test size(result["simulation_result"], 2) == 21 + +# Scientific test. Check that the model in [Kimchi et al. 2020](https://www.science.org/doi/10.1126/sciadv.abc1939) produces oscillatory solutions +# TODO: Postponing until runtime is reduced and simulations validated +""" +connectivity = [1 -1;1 -1] +samples = 5000 +hyperparameters = deepcopy(DEFAULT_SEN_HYPERPARAMETERS) +hyperparameters["random_seed"] = 1234 +hyperparameters["parameter_limits"] = Dict( + "k_b" => [1e-2, 1e0], + "k_u" => [1e-3, 1e3], + "n" => [2.0, 2.0], + "κ_tot" => [1e-3, 1e2], + "η" => [1e-1, 1e1], + "P̃" => [1e-8, 1e-2], +) +hyperparameters["sampling_scales"] = Dict( + "k_b" => "log", + "k_u" => "log", + "n" => "linear", + "κ_tot" => "log", + "η" => "log", + "P̃" => "log", +) +result = find_sen_oscillations(connectivity, samples; + hyperparameters=hyperparameters) +# TODO: Add test +""" \ No newline at end of file From e53e0d4771a9599cc4adeb3db86b38712308eb56 Mon Sep 17 00:00:00 2001 From: ftavella Date: Thu, 25 Jan 2024 13:59:47 -0500 Subject: [PATCH 3/6] Fix issues related to package updates --- src/split_enzyme_network.jl | 2 +- test/split_enzyme_network/sen_test.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/split_enzyme_network.jl b/src/split_enzyme_network.jl index 982aa20..63f5aa6 100644 --- a/src/split_enzyme_network.jl +++ b/src/split_enzyme_network.jl @@ -90,7 +90,7 @@ function split_enzyme_network(connectivity::AbstractMatrix) # Does the network need a constitutive kinase? constitutive_kinase_needed = needs_constitutive_kinase(connectivity) @variables t - @species (X(t))[collect(1:N)] + @species (X(t))[1:N] if constitutive_kinase_needed @parameters k_b[1:N], k_u[1:N], n[1:N], κ_tot[1:N], η, K̃, P̃ else diff --git a/test/split_enzyme_network/sen_test.jl b/test/split_enzyme_network/sen_test.jl index 87f253f..4fb811f 100644 --- a/test/split_enzyme_network/sen_test.jl +++ b/test/split_enzyme_network/sen_test.jl @@ -59,7 +59,7 @@ result = BiologicalOscillations.needs_constitutive_kinase(doesnt_need_kinase_3) @test_throws DomainError split_enzyme_network([0 1;3 -1]) ## Correct reactions - Repressilator @variables t -@species (X(t))[collect(1:3)] +@species (X(t))[1:3] @parameters k_b[1:3], k_u[1:3], n[1:3], κ_tot[1:3], η, P̃ true_reactions = [] From 334834a85d12756bfd2a1f47d8b875513d1131a7 Mon Sep 17 00:00:00 2001 From: ftavella Date: Thu, 25 Jan 2024 18:25:02 -0500 Subject: [PATCH 4/6] Fix file name typo --- test/split_enzyme_network/{sen_test.jl => sen_tests.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/split_enzyme_network/{sen_test.jl => sen_tests.jl} (100%) diff --git a/test/split_enzyme_network/sen_test.jl b/test/split_enzyme_network/sen_tests.jl similarity index 100% rename from test/split_enzyme_network/sen_test.jl rename to test/split_enzyme_network/sen_tests.jl From 0d82fe4d813daf92ac6620d9848d8c9f438f35de Mon Sep 17 00:00:00 2001 From: ftavella Date: Tue, 30 Jan 2024 14:04:53 -0500 Subject: [PATCH 5/6] Fix CI issues --- src/default_hyperparameters.jl | 2 +- test/split_enzyme_network/sen_tests.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/default_hyperparameters.jl b/src/default_hyperparameters.jl index 1d33d46..50ad51d 100644 --- a/src/default_hyperparameters.jl +++ b/src/default_hyperparameters.jl @@ -174,7 +174,7 @@ Default parameter limits for the generation of sen parameter sets. Keys are the """ DEFAULT_SEN_PARAMETER_LIMITS = Dict( "k_b" => (1e-2, 1e0), - "k_u" => (1e-3, 1e-3), + "k_u" => (1e-3, 1e3), "n" => (2.0, 4.0), "κ_tot" => (1e-3, 1e2), "η" => (1e-1, 1e1), diff --git a/test/split_enzyme_network/sen_tests.jl b/test/split_enzyme_network/sen_tests.jl index 4fb811f..f9c8c5b 100644 --- a/test/split_enzyme_network/sen_tests.jl +++ b/test/split_enzyme_network/sen_tests.jl @@ -201,9 +201,9 @@ for (idx, row) in enumerate(eachrow(all_params)) end # Test that parameters are truly oscillatory and non-oscillatory -# Values taken from run on 01/11/2024 -samples = 1500 -random_seed = 9283 +# Values taken from run on 01/29/2024 +samples = 10000 +random_seed = 1832 connectivity_T2_12 = [0 0 1; 1 0 0; 0 -1 1] hyperparameters = deepcopy(DEFAULT_SEN_HYPERPARAMETERS) hyperparameters["random_seed"] = random_seed @@ -230,7 +230,7 @@ for selected_idx in 1:5 peak_number_found[selected_idx] = peak_num end -expected_peaks = [51, 5, 569, 14, 3] +expected_peaks = [1770.0, 185.0, 285.0, 2335.0, 19.0] @test peak_number_found == expected_peaks # Test that parameters are truly non-oscillatory non_oscillatory_params = result["parameter_sets"]["non_oscillatory"] From eccafeadb02fbea6bec3f32e9c30b5af2f424c8c Mon Sep 17 00:00:00 2001 From: ftavella Date: Tue, 30 Jan 2024 15:18:30 -0500 Subject: [PATCH 6/6] Remove nightly from CI --- .github/workflows/CI.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0421f56..d6b2f18 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -19,7 +19,6 @@ jobs: matrix: version: - '1.8' - - 'nightly' os: - ubuntu-latest arch: