Quantification of time scale separation using sojourn time in nullclines

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

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
+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

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
+@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]