Merge pull request #29 from DenisTitovLab/more-edits-to-prevent-param…

…-snooping-in-qssa-mode

add filtering of qssa removal codes at each iterations

src/data_driven_rate_equation_selection.jl

@@ -95,11 +95,13 @@ function data_driven_rate_equation_selection(
             if forward_model_selection
                 nt_param_removal_codes = forward_selection_next_param_removal_codes(
                     nt_previous_param_removal_codes,
+                    metab_names,
                     num_alpha_params,
                 )
             elseif !forward_model_selection
                 nt_param_removal_codes = reverse_selection_next_param_removal_codes(
                     nt_previous_param_removal_codes,
+                    metab_names,
                     num_alpha_params,
                 )
             end
@@ -345,33 +347,11 @@ function calculate_all_parameter_removal_codes_w_num_params(
     end
     nt_param_removal_codes =
         [NamedTuple{param_removal_code_names}(x) for x in unique(codes_with_num_params)]
-
-    # ensure that if K_S1 = Inf then all K_S1_S2 and all other K containing S1 in qssa cannot be 2, which stands for (K_S1_S2)^2 = K_S1 * K_S2
-    if any([occursin("allo", string(key)) for key in keys(nt_param_removal_codes[1])])
-        filtered_nt_param_removal_codes = nt_param_removal_codes
-    else
-        filtered_nt_param_removal_codes = NamedTuple[]
-        for nt_param_removal_code in nt_param_removal_codes
-            if all(
-                nt_param_removal_code[Symbol("K_" * string(metab))] != 1 for
-                metab in metab_names
-            )
-                push!(filtered_nt_param_removal_codes, nt_param_removal_code)
-            else
-                one_metab_codes = metab_names[findall(
-                    nt_param_removal_code[Symbol("K_" * string(metab))] == 1 for
-                    metab in metab_names
-                )]
-                if all(
-                    nt_param_removal_code[param_name] != 2 for
-                    param_name in keys(nt_param_removal_code) if
-                    any(occursin.(string.(one_metab_codes), string(param_name)))
-                )
-                    push!(filtered_nt_param_removal_codes, nt_param_removal_code)
-                end
-            end
-        end
-    end
+    filtered_nt_param_removal_codes =
+        filter_param_removal_codes_to_prevent_wrong_param_combos(
+            nt_param_removal_codes,
+            metab_names,
+        )
     return filtered_nt_param_removal_codes
 end
 
@@ -440,6 +420,7 @@ Calculate `nt_param_removal_codes` with `num_params` including non-zero term com
 """
 function forward_selection_next_param_removal_codes(
     nt_previous_param_removal_codes::Vector{T} where {T<:NamedTuple},
+    metab_names::Tuple{Symbol,Vararg{Symbol}},
     num_alpha_params::Int,
 )
     param_removal_code_names = keys(nt_previous_param_removal_codes[1])
@@ -473,14 +454,20 @@ function forward_selection_next_param_removal_codes(
     end
     nt_param_removal_codes =
         [NamedTuple{param_removal_code_names}(x) for x in unique(next_param_removal_codes)]
-    return nt_param_removal_codes
+    filtered_nt_param_removal_codes =
+        filter_param_removal_codes_to_prevent_wrong_param_combos(
+            nt_param_removal_codes,
+            metab_names,
+        )
+    return filtered_nt_param_removal_codes
 end
 
 """
 Use `nt_previous_param_removal_codes` to calculate `nt_next_param_removal_codes` that have one additional zero elements except for for elements <= `num_alpha_params` from the end
 """
 function reverse_selection_next_param_removal_codes(
     nt_previous_param_removal_codes::Vector{T} where {T<:NamedTuple},
+    metab_names::Tuple{Symbol,Vararg{Symbol}},
     num_alpha_params::Int,
 )
     param_removal_code_names = keys(nt_previous_param_removal_codes[1])
@@ -497,5 +484,44 @@ function reverse_selection_next_param_removal_codes(
     end
     nt_param_removal_codes =
         [NamedTuple{param_removal_code_names}(x) for x in unique(next_param_removal_codes)]
-    return nt_param_removal_codes
+    filtered_nt_param_removal_codes =
+        filter_param_removal_codes_to_prevent_wrong_param_combos(
+            nt_param_removal_codes,
+            metab_names,
+        )
+    return filtered_nt_param_removal_codes
+end
+
+"""Filter removal codes to ensure that if K_S1 = Inf then all K_S1_S2 and all other K containing S1 in qssa cannot be 2, which stands for (K_S1_S2)^2 = K_S1 * K_S2"""
+function filter_param_removal_codes_to_prevent_wrong_param_combos(
+    nt_param_removal_codes,
+    metab_names::Tuple{Symbol,Vararg{Symbol}}
+)
+    # ensure that if K_S1 = Inf then all K_S1_S2 and all other K containing S1 in qssa cannot be 2, which stands for (K_S1_S2)^2 = K_S1 * K_S2
+    if any([occursin("allo", string(key)) for key in keys(nt_param_removal_codes[1])])
+        filtered_nt_param_removal_codes = nt_param_removal_codes
+    else
+        filtered_nt_param_removal_codes = NamedTuple[]
+        for nt_param_removal_code in nt_param_removal_codes
+            if all(
+                nt_param_removal_code[Symbol("K_" * string(metab))] != 1 for
+                metab in metab_names
+            )
+                push!(filtered_nt_param_removal_codes, nt_param_removal_code)
+            else
+                one_metab_codes = metab_names[findall(
+                    nt_param_removal_code[Symbol("K_" * string(metab))] == 1 for
+                    metab in metab_names
+                )]
+                if all(
+                    nt_param_removal_code[param_name] != 2 for
+                    param_name in keys(nt_param_removal_code) if
+                    any(occursin.(string.(one_metab_codes), string(param_name)))
+                )
+                    push!(filtered_nt_param_removal_codes, nt_param_removal_code)
+                end
+            end
+        end
+    end
+    return filtered_nt_param_removal_codes
 end

test/tests_for_optimal_rate_eq_selection.jl

@@ -10,6 +10,7 @@ using BenchmarkTools
 
 #test forward_selection_next_param_removal_codes
 num_metabolites = rand(4:8)
+metab_names = Tuple(Symbol("Metabolite$(i)") for i = 1:num_metabolites)
 n_alphas = rand(1:4)
 num_previous_step_params = rand((2+num_metabolites):(3+2*num_metabolites))
 num_params = num_previous_step_params - 1
@@ -44,6 +45,7 @@ param_removal_code_names
 nt_funct_output_param_subset_codes =
     DataDrivenEnzymeRateEqs.forward_selection_next_param_removal_codes(
         nt_previous_param_removal_codes,
+        metab_names,
         n_alphas,
     )
 funct_output_param_subset_codes = [values(nt) for nt in nt_funct_output_param_subset_codes]
@@ -91,6 +93,7 @@ end
 
 #test reverse_selection_next_param_removal_codes
 num_metabolites = rand(4:8)
+metab_names = Tuple(Symbol("Metabolite$(i)") for i = 1:num_metabolites)
 n_alphas = rand(1:4)
 num_previous_step_params = rand((2+num_metabolites):(3+2*num_metabolites))
 num_params = num_previous_step_params + 1
@@ -123,6 +126,7 @@ nt_previous_param_removal_codes =
 nt_funct_output_param_subset_codes =
     DataDrivenEnzymeRateEqs.reverse_selection_next_param_removal_codes(
         nt_previous_param_removal_codes,
+        metab_names,
         n_alphas,
     )
 funct_output_param_subset_codes = [values(nt) for nt in nt_funct_output_param_subset_codes]
@@ -229,25 +233,59 @@ unidentifiable_params =
         :alpha_ATP_Phenylalanine,
     ) for unidentifiable_param in unidentifiable_params
 )
+
+#test filter_param_removal_codes_to_prevent_wrong_param_combos
+nt_param_removal_codes = [
+    NamedTuple{(:Vmax, :K_S1, :K_S2, :K_S1_S2)}(combo) for
+    combo in Iterators.product([[0, 1], [0, 1], [0, 1], [0, 1, 2]]...)
+]
+metab_names = (:S1, :S2)
+filtered_nt =
+    DataDrivenEnzymeRateEqs.filter_param_removal_codes_to_prevent_wrong_param_combos(
+        nt_param_removal_codes,
+        metab_names,
+    )
+correct_answer = [
+    (Vmax = 0, K_S1 = 0, K_S2 = 0, K_S1_S2 = 0)
+    (Vmax = 1, K_S1 = 0, K_S2 = 0, K_S1_S2 = 0)
+    (Vmax = 0, K_S1 = 1, K_S2 = 0, K_S1_S2 = 0)
+    (Vmax = 1, K_S1 = 1, K_S2 = 0, K_S1_S2 = 0)
+    (Vmax = 0, K_S1 = 0, K_S2 = 1, K_S1_S2 = 0)
+    (Vmax = 1, K_S1 = 0, K_S2 = 1, K_S1_S2 = 0)
+    (Vmax = 0, K_S1 = 1, K_S2 = 1, K_S1_S2 = 0)
+    (Vmax = 1, K_S1 = 1, K_S2 = 1, K_S1_S2 = 0)
+    (Vmax = 0, K_S1 = 0, K_S2 = 0, K_S1_S2 = 1)
+    (Vmax = 1, K_S1 = 0, K_S2 = 0, K_S1_S2 = 1)
+    (Vmax = 0, K_S1 = 1, K_S2 = 0, K_S1_S2 = 1)
+    (Vmax = 1, K_S1 = 1, K_S2 = 0, K_S1_S2 = 1)
+    (Vmax = 0, K_S1 = 0, K_S2 = 1, K_S1_S2 = 1)
+    (Vmax = 1, K_S1 = 0, K_S2 = 1, K_S1_S2 = 1)
+    (Vmax = 0, K_S1 = 1, K_S2 = 1, K_S1_S2 = 1)
+    (Vmax = 1, K_S1 = 1, K_S2 = 1, K_S1_S2 = 1)
+    (Vmax = 0, K_S1 = 0, K_S2 = 0, K_S1_S2 = 2)
+    (Vmax = 1, K_S1 = 0, K_S2 = 0, K_S1_S2 = 2)
+]
+@test filtered_nt == correct_answer
+
+
 #Load and process data
 LDH_data_for_fit = CSV.read(joinpath(@__DIR__, "Data_for_tests/LDH_data.csv"), DataFrame)
 #Add source column that uniquely identifies a figure from publication
 LDH_data_for_fit.source = LDH_data_for_fit.Article .* "_" .* LDH_data_for_fit.Fig
 data = LDH_data_for_fit
 LDH_enzyme = (;
-    substrates=[:Pyruvate, :NADH],
-    products=[:Lactate, :NAD],
-    regulators=[],
-    Keq=20_000.0,
-    rate_equation_name=:ldh_rate_equation,
+    substrates = [:Pyruvate, :NADH],
+    products = [:Lactate, :NAD],
+    regulators = [],
+    Keq = 20_000.0,
+    rate_equation_name = :ldh_rate_equation,
 )
 metab_names, param_names = @derive_general_qssa_rate_eq(LDH_enzyme)
 unidentifiable_params =
     DataDrivenEnzymeRateEqs.find_practically_unidentifiable_params(data, param_names)
 @test all(
-    unidentifiable_param ∈ (
-        :K_Pyruvate_NADH_Lactate_NAD,
-    ) for unidentifiable_param in unidentifiable_params
+    unidentifiable_param ∈ (:K_Pyruvate_NADH_Lactate_NAD,) for
+    unidentifiable_param in unidentifiable_params
 )