diff --git a/src/BiologicalOscillations.jl b/src/BiologicalOscillations.jl index 39fb39d..7e9d7cd 100644 --- a/src/BiologicalOscillations.jl +++ b/src/BiologicalOscillations.jl @@ -10,7 +10,7 @@ export elowitz_2000, guan_2008 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 +export calculate_main_frequency, calculate_amplitude, is_ODE_oscillatory, calculate_sojourn_time_fractions_in_nullclines # Protein interaction network export protein_interaction_network, pin_parameters, pin_timescale, pin_parameter_sets export pin_equilibration_times, find_pin_oscillations, pin_nodes_edges diff --git a/src/feature_calculation.jl b/src/feature_calculation.jl index e5df206..c6702eb 100644 --- a/src/feature_calculation.jl +++ b/src/feature_calculation.jl @@ -135,4 +135,75 @@ function is_ODE_oscillatory(frequency_data::Dict, amplitude_data::Dict; freq_var else return false end -end \ No newline at end of file +end + + +""" + calculate_sojourn_time_fractions_in_nullclines(ode_solution::ODESolution) + +Calculate the time fractions that the dynamics stays near nullclines. These numbers can be use as a metric to quantify the time scale separation of oscillations. + +# Arguments (Required) +- `ode_solution::ODESolution`: The full output of the `solve()` function from `DifferentialEquations.jl`. The solution should be oscillatory and equilibrated at least for the last cycle. The solution should also be long enough for accurate determination of periods. + +# Arguments (Optional) +- `eps::Float64`: Threshold to determine if trajectories are close enough to nullclines. + +# Returns +- `tss::Vector{Float64}`: Each element is the sojourn time fraction (sojourn time per period / period) computed for the respective nullcline. The length of the output vector equals the ODE dimension. +""" +function calculate_sojourn_time_fractions_in_nullclines(ode_solution::ODESolution; eps::Float64=1e-2) + period = 1 / mean(calculate_main_frequency(ode_solution, length(ode_solution.t), length(ode_solution.t))["frequency"]) + + idx_last_cycle = ode_solution.t .> ode_solution.t[end] - period + idx_last_cycle[findfirst(diff(idx_last_cycle) .== 1)] = 1 + # add one more data point to fully complete a cycle + + function gradient(v, grid) + # 1D equivalent of numpy.gradient with default inputs + # (second order accurate finite difference for interior and first order for edges) + grad = zeros(length(v)) + + d = diff(grid) + + h_s = d[1:end - 1] + h_d = d[2:end] + + f_plus = v[3:end] + f_0 = v[2:end - 1] + f_minus = v[1:end - 2] + + # second order central difference for interior points + grad[2:end - 1] = @. (h_s ^ 2 * f_plus + (h_d ^ 2 - h_s ^ 2) * f_0 - h_d ^ 2 * f_minus) / (h_s * h_d * (h_d + h_s)) + + # first order forward/backward differences for boundary points + grad[1] = (v[2] - v[1]) / d[1] + grad[end] = (v[end] - v[end - 1]) / d[end] + + return grad + end + + N = length(ode_solution.u[1]) + tss = [] + + for i = 1:N + t = ode_solution.t[idx_last_cycle] + sol = [v[i] for v in ode_solution.u[idx_last_cycle]] + + sol = (sol .- minimum(sol)) / (maximum(sol) - minimum(sol)) + + grad = gradient(sol, t) + + # thresholding to determine how close trajectories are to slow manifolds + nc = abs.(grad) .< maximum(abs.(grad)) * eps + + nc = (nc[1:end - 1] + nc[2:end]) / 2 + + tss_ = sum(nc .* diff(t)) / (maximum(t) - minimum(t)) + # total time that the trajectory stays near i-th variable nullcline + + push!(tss, tss_) + end + + return tss +end diff --git a/test/feature_calculation_tests.jl b/test/feature_calculation_tests.jl index 80457d1..9cc9376 100644 --- a/test/feature_calculation_tests.jl +++ b/test/feature_calculation_tests.jl @@ -139,4 +139,121 @@ frequency_data = calculate_main_frequency(repressilator_sol, signal_sampling, sp @test mean(frequency_data["power"]) < 1e-5 @test all(isnan.(amplitude_data["amplitude"])) -@test is_ODE_oscillatory(frequency_data, amplitude_data) == false \ No newline at end of file +@test is_ODE_oscillatory(frequency_data, amplitude_data) == false + +# time scale separation + +# sojourn time in nullclines, calculate_sojourn_time_fractions_in_nullclines +# parameter set #3 and #4 give spiky solutions + +connectivity_interaction = [ + [1 0 -1; 1 0 0; 0 1 0], + [1 0 -1; 1 0 0; 0 1 0], + [1 0 -1; 1 0 0; 0 1 0], + [1 0 -1; 1 0 0; 0 1 0], + [1 0 -1; 1 0 0; 0 1 0], + [0 0 -1; 1 0 0; 0 1 0], + [0 0 -1; 1 0 0; 0 1 0], + [0 0 -1; 1 0 0; 0 1 0], + [0 0 -1; 1 0 0; 0 1 0], + [0 0 -1; 1 0 0; 0 1 0], +] + +α = [ + [1,0.215443,0.284804], + [1,0.0278256,5.09414], + [1,0.0159228,2.65609], + [1,0.033516,0.453488], + [1,0.599484,0.010975], + [1,1.21958,2.1234], + [1,0.0494845,0.0246626], + [1,0.0214404,0.0779272], + [1,0.135554,0.260941], + [1,1.3571,0.0335469], +] + +β = [ + [0.0305386,32.7455,3.51119], + [0.070548,24.7708,17.0735], + [14.1747,35.9381,4.64159], + [29.8365,15.5568,57.2237], + [1.14976,10.7227,8.90215], + [0.0214009,11.8116,32.565], + [0.06375,14.8155,30.3353], + [0.614015,6.77147,24.8164], + [6.09647,19.4057,37.1156], + [0.0578197,30.6725,5.32933], +] + +γ = [ + [159.228,291.505,7220.81,642.807], + [559.081,422.924,1555.68,2364.49], + [4977.02,2983.65,739.072,231.013], + [1629.75,8302.18,4328.76,126.186], + [486.26,9545.48,1291.55,774.264], + [396.698,5055.57,234.119], + [672.485,215.791,1255.19], + [5447.2,8038.79,1207.56], + [942.969,5462.28,1565.74], + [336.863,4114.92,161.221], +] + +κ = [ + [0.975758,0.450505,0.531313,0.765657], + [0.862626,0.951515,0.692929,0.79798], + [0.660606,0.733333,0.644444,0.749495], + [0.668687,0.256566,0.507071,0.264646], + [0.59596,0.216162,0.353535,0.757576], + [0.824862,0.360336,0.816622], + [0.526673,0.90015,0.586039], + [0.827583,0.534513,0.589239], + [0.631803,0.69861,0.20464], + [0.947435,0.246405,0.524912], +] + +η = [ + [4.79798,3.62626,3.90909,3.26263], + [3.0202,4.75758,3.74747,3.54545], + [1.88889,4.31313,1.28283,3.70707], + [2.33333,4.07071,3.42424,3.74747], + [4.11111,4.47475,2.29293,4.39394], + [4.33993,3.81748,2.94179], + [2.13491,4.39114,4.79118], + [3.78468,3.38704,2.15292], + [4.46355,2.20132,4.17912], + [2.50375,3.12861,4.76478], +] + +tspan = (0.0, 2.0) +initial_condition = [0.6, 0.6, 0.6] + +max_sojourn_time = [] + +for (i, c) in enumerate(connectivity_interaction) + model_interaction = protein_interaction_network(c) + params = pin_parameters(model_interaction, α[i], β[i], γ[i], κ[i], η[i]) + equilibration = ODEProblem(model_interaction, initial_condition, tspan, params) + sol_eq = solve(equilibration) + ode_problem = ODEProblem(model_interaction, sol_eq.u[end], tspan, params) + sol = solve(ode_problem) + + push!(max_sojourn_time, maximum(calculate_sojourn_time_fractions_in_nullclines(sol))) +end + +@test max_sojourn_time[3] > max_sojourn_time[1] +@test max_sojourn_time[3] > max_sojourn_time[2] +@test max_sojourn_time[3] > max_sojourn_time[5] +@test max_sojourn_time[3] > max_sojourn_time[6] +@test max_sojourn_time[3] > max_sojourn_time[7] +@test max_sojourn_time[3] > max_sojourn_time[8] +@test max_sojourn_time[3] > max_sojourn_time[9] +@test max_sojourn_time[3] > max_sojourn_time[10] + +@test max_sojourn_time[4] > max_sojourn_time[1] +@test max_sojourn_time[4] > max_sojourn_time[2] +@test max_sojourn_time[4] > max_sojourn_time[5] +@test max_sojourn_time[4] > max_sojourn_time[6] +@test max_sojourn_time[4] > max_sojourn_time[7] +@test max_sojourn_time[4] > max_sojourn_time[8] +@test max_sojourn_time[4] > max_sojourn_time[9] +@test max_sojourn_time[4] > max_sojourn_time[10]