From 72a1c585c4f2542a461a97cccf25544d3b53c4d3 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Mon, 28 Oct 2024 13:27:28 +0000 Subject: [PATCH] build based on 8fcc44c --- dev/.documenter-siteinfo.json | 2 +- dev/1-parametric-models.jl | 171 +++++++------- dev/1-parametric-models/index.html | 35 +-- dev/2-differentiate-enzyme-model.jl | 58 +++-- dev/2-differentiate-enzyme-model/index.html | 238 ++++++++++---------- dev/3-parameter-estimation.jl | 34 +-- dev/3-parameter-estimation/98c59523.png | Bin 0 -> 32477 bytes dev/3-parameter-estimation/ffa91456.png | Bin 34990 -> 0 bytes dev/3-parameter-estimation/index.html | 43 ++-- dev/index.html | 2 +- dev/objects.inv | Bin 1222 -> 1260 bytes dev/reference/index.html | 52 ++--- dev/search_index.js | 2 +- 13 files changed, 319 insertions(+), 318 deletions(-) create mode 100644 dev/3-parameter-estimation/98c59523.png delete mode 100644 dev/3-parameter-estimation/ffa91456.png diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index ad6e78c..5519ea6 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-10-24T17:31:43","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-10-28T13:27:25","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/1-parametric-models.jl b/dev/1-parametric-models.jl index bb7160a..c8007ce 100644 --- a/dev/1-parametric-models.jl +++ b/dev/1-parametric-models.jl @@ -16,7 +16,7 @@ using DifferentiableMetabolism -using Symbolics +using FastDifferentiation using ConstraintTrees using COBREXA using Tulip @@ -46,7 +46,7 @@ base_model.fluxes # ## Add parameters to the model # Make bound of r2 and mass balance of m3 parameters -Symbolics.@variables r2bound m3bound +@variables r2bound m3bound m.fluxes.r2 = ConstraintTrees.Constraint(m.fluxes.r2.value, -2 * ParameterBetween(r2bound, 0)) @@ -55,7 +55,7 @@ m.flux_stoichiometry.m3 = ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2) # # add parametric constraints -Symbolics.@variables p[1:4] +p = make_variables(:p, 4) m *= :linparam^ConstraintTrees.Constraint( @@ -65,11 +65,11 @@ m *= # substitute params into model parameter_substitutions = Dict( - r2bound => 4.0, - m3bound => 0.1, # lose some mass here - p[1] => 1.0, - p[2] => 1.0, - p[3] => 4.0, + :r2bound => 4.0, + :m3bound => 0.1, # lose some mass here + :p1 => 1.0, + :p2 => 1.0, + :p3 => 4.0, ) m_noparams, _, _, _ = optimized_constraints_with_parameters( @@ -85,7 +85,7 @@ m_noparams.fluxes # ## Change the parameters and re-solve # substitute params into model -parameter_substitutions[m3bound] = 0.0 +parameter_substitutions[:m3bound] = 0.0 m_noparams, _, _, _ = optimized_constraints_with_parameters( m, @@ -99,7 +99,7 @@ m_noparams.fluxes # ## Quadratic parameters also work -Symbolics.@variables q[1:6] +q = make_variables(:q, 6) m.objective = ConstraintTrees.Constraint( value = sum( @@ -110,7 +110,8 @@ m.objective = ConstraintTrees.Constraint( m *= :objective_bound^ConstraintTrees.Constraint(value = m.fluxes.r6.value, bound = 2.0) -parameter_substitutions = merge(parameter_substitutions, Dict(zip(q, fill(1.0, 6)))) +parameter_substitutions = + merge(parameter_substitutions, Dict(v.node_value => 1.0 for v in q)) m_noparams, _, _, _ = optimized_constraints_with_parameters( m, @@ -121,138 +122,136 @@ m_noparams, _, _, _ = optimized_constraints_with_parameters( ) m_noparams.fluxes +v(x) = DifferentiableMetabolism.substitute(x, _ -> throw("no subst!")) #src @test isapprox(m_noparams.objective, 11.0; atol = TEST_TOLERANCE) #src pb = ParameterBetween(1, 2.0) #src -@test (-pb).lower == -1 && (-pb).upper == -2 #src -@test (pb * 2).lower == 2 && (pb * 2).upper == 4 #src -@test (pb / 2).lower == 0.5 && (pb / 2).upper == 1 #src +@test v((-pb).lower) == -1 && v((-pb).upper) == -2 #src +@test v((pb * 2).lower) == 2 && v((pb * 2).upper) == 4 #src +@test v((pb / 2).lower) == 0.5 && v((pb / 2).upper) == 1 #src pe = ParameterEqualTo(1) #src -@test (-pe).equal_to == -1 #src -@test (pe * 2).equal_to == 2 #src -@test (pe / 2).equal_to == 0.5 #src +@test v((-pe).equal_to) == -1 #src +@test v((pe * 2).equal_to) == 2 #src +@test v((pe / 2).equal_to) == 0.5 #src plv1 = ParameterLinearValue(ConstraintTrees.LinearValue(1)) #src -@test all(plv1.idxs .== [0]) && all(plv1.weights .== [Num(1.0)]) #src +@test plv1.idxs == [0] && all(v.(plv1.weights) .== 1.0) #src plv2 = ParameterLinearValue([0], [1.0]) #src -@test all(plv2.idxs .== [0]) && all(plv2.weights .== [Num(1.0)]) #src +@test plv2.idxs == [0] && all(v.(plv2.weights) .== 1.0) #src plv3 = ParameterLinearValue(0) #src @test isempty(plv3.idxs) && isempty(plv3.weights) #src plv4 = ParameterLinearValue(1) #src -@test all(plv4.idxs .== [0]) && all(plv4.weights .== [Num(1.0)]) #src +@test plv4.idxs == [0] && all(v.(plv4.weights) .== 1.0) #src plv5 = convert(ParameterLinearValue, 1) #src -@test all(plv5.idxs .== [0]) && all(plv5.weights .== [Num(1.0)]) #src +@test plv5.idxs == [0] && all(v.(plv5.weights) .== 1.0) #src plv6 = zero(ParameterLinearValue) #src @test isempty(plv6.idxs) && isempty(plv6.weights) #src plv7 = plv5 + 1 #src -@test all(plv7.idxs .== [0]) && all(plv7.weights .== [Num(2.0)]) #src +@test all(v.(plv7.idxs) .== 0) && all(v.(plv7.weights) .== 2.0) #src plv8 = ParameterLinearValue([1, 2, 3], [1, 2, 3]) #src plv9 = ParameterLinearValue([1, 3], [4.0, 5.0]) #src plv10 = plv8 + plv9 #src -@test all(plv10.idxs .== [1, 2, 3]) && all(plv10.weights .== [Num(5), Num(2), Num(8)]) #src +@test plv10.idxs == [1, 2, 3] && v.(plv10.weights) == [5, 2, 8] #src plv10 = plv9 + plv8 #src -@test all(plv10.idxs .== [1, 2, 3]) && all(plv10.weights .== [Num(5), Num(2), Num(8)]) #src +@test plv10.idxs == [1, 2, 3] && v.(plv10.weights) == [5, 2, 8] #src plv11 = plv5 - 2 #src -@test all(plv11.idxs .== [0]) && all(plv11.weights .== [Num(-1)]) #src +@test plv11.idxs == [0] && v.(plv11.weights) == [-1] #src plv12 = 2 - plv5 #src -@test all(plv12.idxs .== [0]) && all(plv12.weights .== [Num(1)]) #src +@test plv12.idxs == [0] && v.(plv12.weights) == [1] #src plv13 = plv5 / -2 #src -@test all(plv13.idxs .== [0]) && all(plv13.weights .== [Num(-0.5)]) #src +@test plv13.idxs == [0] && v.(plv13.weights) == [-0.5] #src plv14 = plv8 - plv9 #src -@test all(plv14.idxs .== [1, 2, 3]) && all(plv14.weights .== [Num(-3), Num(2), Num(-2)]) #src +@test plv14.idxs == [1, 2, 3] && v.(plv14.weights) == [-3, 2, -2] #src plv15 = 1.0 + ParameterLinearValue(0) #src -@test all(plv15.idxs .== [0]) && all(plv15.weights .== [Num(1)]) #src +@test plv15.idxs == [0] && v.(plv15.weights) == [1] #src pqv1 = ParameterQuadraticValue([(1, 1)], [1]) #src pqv2 = ParameterQuadraticValue(ConstraintTrees.QuadraticValue([(1, 1)], [1])) #src -@test all(pqv2.idxs[1] .== (1, 1)) && all(pqv2.weights .== Num(1)) #src +@test pqv2.idxs == [(1, 1)] && v.(pqv2.weights) == [1] #src pqv3 = ParameterQuadraticValue(0) #src @test isempty(pqv3.idxs) && isempty(pqv3.weights) #src pqv4 = ParameterQuadraticValue(ParameterLinearValue([1], [1])) #src -@test all(pqv4.idxs .== [(0, 1)]) && all(pqv4.weights .== [Num(1)]) #src +@test pqv4.idxs == [(0, 1)] && v.(pqv4.weights) == [1] #src pqv5 = ParameterQuadraticValue(1) #src -@test all(pqv5.idxs .== [(0, 0)]) && all(pqv5.weights .== [Num(1)]) #src +@test pqv5.idxs == [(0, 0)] && v.(pqv5.weights) == [1] #src pqv6 = convert(ParameterQuadraticValue, 1) #src -@test all(pqv6.idxs .== [(0, 0)]) && all(pqv6.weights .== [Num(1)]) #src +@test pqv6.idxs == [(0, 0)] && v.(pqv6.weights) == [1] #src pqv7 = convert(ParameterQuadraticValue, ParameterLinearValue([1], [1])) #src -@test all(pqv7.idxs .== [(0, 1)]) && all(pqv7.weights .== [Num(1)]) #src +@test pqv7.idxs == [(0, 1)] && v.(pqv7.weights) == [1] #src pqv8 = zero(ParameterQuadraticValue) #src @test isempty(pqv8.idxs) && isempty(pqv8.weights) #src pqv9 = 1 + pqv1 #src -@test all(pqv9.idxs .== [(0, 0), (1, 1)]) && all(pqv9.weights .== [Num(1), Num(1)]) #src +@test pqv9.idxs == [(0, 0), (1, 1)] && v.(pqv9.weights) == [1, 1] #src pqv10 = ParameterLinearValue([1], [1]) + pqv1 #src -@test all(pqv10.idxs .== [(0, 1), (1, 1)]) && all(pqv10.weights .== [Num(1), Num(1)]) #src +@test pqv10.idxs == [(0, 1), (1, 1)] && v.(pqv10.weights) == [1, 1] #src pqv11 = 3 - pqv1 #src -@test all(pqv11.idxs .== [(0, 0), (1, 1)]) && all(pqv11.weights .== [Num(3), Num(-1)]) #src +@test pqv11.idxs == [(0, 0), (1, 1)] && v.(pqv11.weights) == [3, -1] #src pqv12 = ParameterLinearValue([1], [1]) - pqv1 #src -@test all(pqv12.idxs .== [(0, 1), (1, 1)]) && all(pqv12.weights .== [Num(1), Num(-1)]) #src +@test pqv12.idxs == [(0, 1), (1, 1)] && v.(pqv12.weights) == [1, -1] #src pqv13 = pqv1 - ParameterLinearValue([1], [1]) #src -@test all(pqv13.idxs .== [(0, 1), (1, 1)]) && all(pqv13.weights .== [Num(-1), Num(1)]) #src +@test pqv13.idxs == [(0, 1), (1, 1)] && v.(pqv13.weights) == [-1, 1] #src pqv14 = 42 * pqv1 #src -@test all(pqv14.idxs .== [(1, 1)]) && all(pqv14.weights .== [Num(42)]) #src +@test pqv14.idxs == [(1, 1)] && v.(pqv14.weights) == [42] #src pqv15 = pqv1 - (3 * pqv1) #src -@test all(pqv15.idxs .== [(1, 1)]) && all(pqv15.weights .== [Num(-2)]) #src +@test pqv15.idxs == [(1, 1)] && v.(pqv15.weights) == [-2] #src pqv16 = pqv1 / 2 #src -@test all(pqv16.idxs .== [(1, 1)]) && all(pqv16.weights .== [Num(0.5)]) #src -pqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src -@test all(pqv17.idxs .== [(1, 2), (2, 2)]) && all(pqv17.weights .== [Num(2), Num(1)]) #src +@test pqv16.idxs == [(1, 1)] && v.(pqv16.weights) == [0.5] #src +pqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src +@test pqv17.idxs == [(1, 2), (2, 2)] && v.(pqv17.weights) == [2, 1] #src pqv18 = ParameterQuadraticValue([(1, 1), (2, 2)], [1.0, 3.0]) #src pqv19 = pqv18 + pqv1 #src -@test all(pqv19.idxs .== [(1, 1), (2, 2)]) && all(pqv19.weights .== [Num(2), Num(3)]) #src +@test pqv19.idxs == [(1, 1), (2, 2)] && v.(pqv19.weights) == [2, 3] #src pqv20 = pqv1 + pqv18 #src -@test all(pqv20.idxs .== [(1, 1), (2, 2)]) && all(pqv19.weights .== [Num(2), Num(3)]) #src +@test pqv20.idxs == [(1, 1), (2, 2)] && v.(pqv19.weights) == [2, 3] #src pqv21 = ParameterQuadraticValue([(1, 1), (2, 2)], [1, 2]) #src -@test all(pqv21.idxs .== [(1, 1), (2, 2)]) && all(pqv21.weights .== [Num(1), Num(2)]) #src +@test pqv21.idxs == [(1, 1), (2, 2)] && v.(pqv21.weights) == [1, 2] #src pqv22 = convert(ParameterQuadraticValue, ConstraintTrees.LinearValue([1], [1])) #src -@test all(pqv22.idxs .== [(0, 1)]) && all(pqv22.weights .== [Num(1)]) #src +@test pqv22.idxs == [(0, 1)] && v.(pqv22.weights) == [1] #src pqv23 = convert( #src ParameterQuadraticValue, #src ConstraintTrees.QuadraticValue([(1, 1), (2, 2)], [1, 2]), #src ) #src -@test all(pqv23.idxs .== [(1, 1), (2, 2)]) && all(pqv23.weights .== [Num(1), Num(2)]) #src -pqv24 = Num(5) - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src -@test all(pqv24.idxs .== [(0, 0), (1, 1)]) && all(pqv24.weights .== [Num(5), Num(-1)]) #src +@test pqv23.idxs == [(1, 1), (2, 2)] && v.(pqv23.weights) == [1, 2] #src +pqv24 = 5 - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src +@test pqv24.idxs == [(0, 0), (1, 1)] && v.(pqv24.weights) == [5, -1] #src pqv25 = pqv18 - 5.0 #src -@test all(pqv25.idxs .== [(0, 0), (1, 1), (2, 2)]) && #src - all(pqv25.weights .== [Num(-5), Num(1), Num(3)]) #src -pr1 = Num(1) + ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr1.idxs .== [0, 1, 2]) && all(pr1.weights .== [Num(1), Num(1), Num(2)]) #src -pr2 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - Num(3) #src -@test all(pr2.idxs .== [0, 1, 2]) && all(pr2.weights .== [Num(-3), Num(1), Num(2)]) #src -pr2b = Num(3) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr2b.idxs .== [0, 1, 2]) && all(pr2b.weights .== [Num(3), Num(-1), Num(-2)]) #src -pr3 = Num(3) * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr3.idxs .== [1, 2]) && all(pr3.weights .== [Num(3), Num(6)]) #src -pr3 = ConstraintTrees.LinearValue([1, 2], [1, 2]) / Num(2) #src -@test all(pr3.idxs .== [1, 2]) && all(pr3.weights .== [Num(0.5), Num(1)]) #src +@test pqv25.idxs == [(0, 0), (1, 1), (2, 2)] && v.(pqv25.weights) == [-5, 1, 3] #src +pr1 = 1 + ConstraintTrees.LinearValue([1, 2], [1, 2]) #src +@test pr1.idxs == [0, 1, 2] && v.(pr1.weights) == [1, 1, 2] #src +pr2 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - 3 #src +@test pr2.idxs == [0, 1, 2] && v.(pr2.weights) == [-3, 1, 2] #src +pr2b = 3 - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src +@test pr2b.idxs == [0, 1, 2] && v.(pr2b.weights) == [3, -1, -2] #src +pr3 = 3 * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src +@test pr3.idxs == [1, 2] && v.(pr3.weights) == [3, 6] #src +pr3 = ConstraintTrees.LinearValue([1, 2], [1, 2]) / 2 #src +@test pr3.idxs == [1, 2] && v.(pr3.weights) == [0.5, 1] #src pr4 = ConstraintTrees.LinearValue([1, 2], [1, 2]) + ParameterLinearValue([1, 2], [1, 2]) #src -@test all(pr4.idxs .== [1, 2]) && all(pr4.weights .== [Num(2), Num(4)]) #src +@test pr4.idxs == [1, 2] && v.(pr4.weights) == [2, 4] #src pr5 = ParameterLinearValue([1, 2], [2, 3]) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr5.idxs .== [1, 2]) && all(pr5.weights .== [Num(1), Num(1)]) #src -pr6 = Num(1) + ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src -@test all(pr6.idxs .== [(0, 0), (1, 1)]) && all(pr6.weights .== [Num(1), Num(1)]) #src -pr7 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) - Num(1) #src -@test all(pr7.idxs .== [(0, 0), (1, 1)]) && all(pr7.weights .== [Num(-1), Num(1)]) #src -pr8 = Num(2) * ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src -@test all(pr8.idxs .== [(1, 1)]) && all(pr8.weights .== [Num(2)]) #src -pr9 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) / Num(2) #src -@test all(pr9.idxs .== [(1, 1)]) && all(pr9.weights .== [Num(0.5)]) #src +@test pr5.idxs == [1, 2] && v.(pr5.weights) == [1, 1] #src +pr6 = 1 + ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src +@test pr6.idxs == [(0, 0), (1, 1)] && v.(pr6.weights) == [1, 1] #src +pr7 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) - 1 #src +@test pr7.idxs == [(0, 0), (1, 1)] && v.(pr7.weights) == [-1, 1] #src +pr8 = 2 * ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src +@test pr8.idxs == [(1, 1)] && v.(pr8.weights) == [2] #src +pr9 = ConstraintTrees.QuadraticValue([(1, 1)], [1]) / 2 #src +@test pr9.idxs == [(1, 1)] && v.(pr9.weights) == [0.5] #src pr10 = #src ConstraintTrees.QuadraticValue([(1, 1)], [1]) + ParameterQuadraticValue([(1, 1)], [1]) #src -@test all(pr10.idxs .== [(1, 1)]) && all(pr10.weights .== [Num(2)]) #src +@test pr10.idxs == [(1, 1)] && v.(pr10.weights) == [2] #src pr11 = #src ParameterQuadraticValue([(1, 1)], [2]) - ConstraintTrees.QuadraticValue([(1, 1)], [1]) #src -@test all(pr11.idxs .== [(1, 1)]) && all(pr11.weights .== [Num(1)]) #src +@test pr11.idxs == [(1, 1)] && v.(pr11.weights) == [1] #src pr12 = ConstraintTrees.LinearValue([1, 2], [1, 2]) + ParameterQuadraticValue([(1, 1)], [2]) #src -@test all(pr12.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src - all(pr12.weights .== [Num(1), Num(2), Num(2)]) #src +@test pr12.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr12.weights) == [1, 2, 2] #src pr13 = ParameterQuadraticValue([(1, 1)], [2]) - ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr13.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src - all(pr13.weights .== [Num(-1), Num(2), Num(-2)]) #src +@test pr13.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr13.weights) == [-1, 2, -2] #src pr14 = ParameterLinearValue([1, 2], [2, 3]) * ConstraintTrees.LinearValue([1, 2], [1, 2]) #src -@test all(pr14.idxs .== [(1, 1), (1, 2), (2, 2)]) && #src - all(pr14.weights .== [Num(2), Num(7), Num(6)]) #src +@test pr14.idxs == [(1, 1), (1, 2), (2, 2)] && v.(pr14.weights) == [2, 7, 6] #src pr15 = ConstraintTrees.LinearValue([1, 2], [1, 2]) - ParameterQuadraticValue([(1, 1)], [2]) #src -@test all(pr15.idxs .== [(0, 1), (1, 1), (0, 2)]) && #src - all(pr15.weights .== [Num(1), Num(-2), Num(2)]) #src -@test ConstraintTrees.substitute( #src - ParameterQuadraticValue([(0, 0), (1, 1)], [1, 2]), #src - [Num(1), Num(2)], #src +@test pr15.idxs == [(0, 1), (1, 1), (0, 2)] && v.(pr15.weights) == [1, -2, 2] #src +@test v( #src + ConstraintTrees.substitute( #src + ParameterQuadraticValue([(0, 0), (1, 1)], [1, 2]), #src + FastDifferentiation.Node.([1.0, 2.0]), #src + ), #src ) == 3 #src diff --git a/dev/1-parametric-models/index.html b/dev/1-parametric-models/index.html index 819926c..7a515e9 100644 --- a/dev/1-parametric-models/index.html +++ b/dev/1-parametric-models/index.html @@ -1,7 +1,7 @@ Parametric constraint-based metabolic models · DifferentiableMetabolism.jl
GitHub

Parametric constraint-based metabolic models

using DifferentiableMetabolism
 
-using Symbolics
+using FastDifferentiation
 using ConstraintTrees
 using COBREXA
 using Tulip
@@ -23,13 +23,13 @@
   :r3 => 1.97554
   :r4 => 1.97554
   :r5 => -0.97554
-  :r6 => 1.0

Add parameters to the model

Make bound of r2 and mass balance of m3 parameters

Symbolics.@variables r2bound m3bound
+  :r6 => 1.0

Add parameters to the model

Make bound of r2 and mass balance of m3 parameters

@variables r2bound m3bound
 
 m.fluxes.r2 =
     ConstraintTrees.Constraint(m.fluxes.r2.value, -2 * ParameterBetween(r2bound, 0))
 
 m.flux_stoichiometry.m3 =
-    ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2)
ConstraintTrees.Constraint(ConstraintTrees.LinearValue([3, 4], [1.0, -1.0]), ParameterEqualTo(0.5m3bound))

add parametric constraints

Symbolics.@variables p[1:4]
+    ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2)
ConstraintTrees.Constraint(ConstraintTrees.LinearValue([3, 4], [1.0, -1.0]), ParameterEqualTo((0.5 * m3bound)))

add parametric constraints

p = make_variables(:p, 4)
 
 m *=
     :linparam^ConstraintTrees.Constraint(
@@ -39,13 +39,13 @@
   :coupling           => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 0 elements =#)
   :flux_stoichiometry => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 4 elements =#)
   :fluxes             => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 6 elements =#)
-  :linparam           => ConstraintTrees.Constraint(ParameterLinearValue(#= ... =#), ParameterBetween(-p[3], 0))
+  :linparam           => ConstraintTrees.Constraint(ParameterLinearValue(#= ... =#), ParameterBetween(-(p3), 0))
   :objective          => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#))

substitute params into model

parameter_substitutions = Dict(
-    r2bound => 4.0,
-    m3bound => 0.1, # lose some mass here
-    p[1] => 1.0,
-    p[2] => 1.0,
-    p[3] => 4.0,
+    :r2bound => 4.0,
+    :m3bound => 0.1, # lose some mass here
+    :p1 => 1.0,
+    :p2 => 1.0,
+    :p3 => 4.0,
 )
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
@@ -60,7 +60,7 @@
   :r3 => 0.05
   :r4 => 4.16342e-11
   :r5 => 3.9
-  :r6 => 3.9

Change the parameters and re-solve

substitute params into model

parameter_substitutions[m3bound] = 0.0
+  :r6 => 3.9

Change the parameters and re-solve

substitute params into model

parameter_substitutions[:m3bound] = 0.0
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
     m,
@@ -74,7 +74,7 @@
   :r3 => 6.43725e-10
   :r4 => 6.45009e-10
   :r5 => 4.0
-  :r6 => 4.0

Quadratic parameters also work

Symbolics.@variables q[1:6]
+  :r6 => 4.0

Quadratic parameters also work

q = make_variables(:q, 6)
 
 m.objective = ConstraintTrees.Constraint(
     value = sum(
@@ -85,7 +85,8 @@
 
 m *= :objective_bound^ConstraintTrees.Constraint(value = m.fluxes.r6.value, bound = 2.0)
 
-parameter_substitutions = merge(parameter_substitutions, Dict(zip(q, fill(1.0, 6))))
+parameter_substitutions =
+    merge(parameter_substitutions, Dict(v.node_value => 1.0 for v in q))
 
 m_noparams, _, _, _ = optimized_constraints_with_parameters(
     m,
@@ -94,6 +95,10 @@
     optimizer = Clarabel.Optimizer,
     sense = Minimal,
 )
-m_noparams.fluxes
-
-pqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src
ParameterQuadraticValue([(1, 2), (2, 2)], Symbolics.Num[2, 1])

This page was generated using Literate.jl.

+m_noparams.fluxes
ConstraintTrees.Tree{Float64} with 6 elements:
+  :r1 => -0.5
+  :r2 => -2.0
+  :r3 => 0.5
+  :r4 => 0.5
+  :r5 => 1.5
+  :r6 => 2.0

This page was generated using Literate.jl.

diff --git a/dev/2-differentiate-enzyme-model.jl b/dev/2-differentiate-enzyme-model.jl index d53e2ea..c2ed938 100644 --- a/dev/2-differentiate-enzyme-model.jl +++ b/dev/2-differentiate-enzyme-model.jl @@ -16,13 +16,13 @@ using DifferentiableMetabolism using AbstractFBCModels -using Symbolics using ConstraintTrees using COBREXA using Tulip using JSONFBCModels import Downloads: download using CairoMakie +using FastDifferentiation !isfile("e_coli_core.json") && download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") @@ -38,12 +38,10 @@ model.reactions[glc_idx]["lower_bound"] = -1000.0 pfl_idx = first(indexin(["PFL"], rids)) model.reactions[pfl_idx]["upper_bound"] = 0.0 -kcats = vcat([Symbolics.@variables $x for x in Symbol.(keys(ecoli_core_reaction_kcats))]...) +rid_kcat = Dict(k => FastDifferentiation.Node(Symbol(k)) for (k,_) in ecoli_core_reaction_kcats) +kcats = Symbol.(keys(ecoli_core_reaction_kcats)) -rid_kcat = Dict(zip(keys(ecoli_core_reaction_kcats), kcats)) - -parameter_values = - Dict(kid => ecoli_core_reaction_kcats[rid] * 3.6 for (rid, kid) in rid_kcat) # change unit to k/h +parameter_values = Dict{Symbol, Float64}() reaction_isozymes = Dict{String,Dict{String,ParameterIsozyme}}() # a mapping from reaction IDs to isozyme IDs to isozyme structs. float_reaction_isozymes = Dict{String,Dict{String,COBREXA.Isozyme}}() #src @@ -52,25 +50,29 @@ for rid in AbstractFBCModels.reactions(model) isnothing(grrs) && continue # skip if no grr available haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available for (i, grr) in enumerate(grrs) + + kcat = ecoli_core_reaction_kcats[rid] * 3.6 # change unit to k/h + parameter_values[Symbol(rid)] = kcat + d = get!(reaction_isozymes, rid, Dict{String,ParameterIsozyme}()) - d2 = get!(float_reaction_isozymes, rid, Dict{String,COBREXA.Isozyme}()) #src - d["isozyme_"*string(i)] = ParameterIsozyme( + d["isozyme_$i"] = ParameterIsozyme( gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes kcat_forward = rid_kcat[rid], kcat_reverse = rid_kcat[rid], ) - d2["isozyme_"*string(i)] = COBREXA.Isozyme( #src + d2 = get!(float_reaction_isozymes, rid, Dict{String,COBREXA.Isozyme}()) #src + d2["isozyme_$i"] = COBREXA.Isozyme( #src gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), #src - kcat_forward = ecoli_core_reaction_kcats[rid] * 3.6, #src - kcat_reverse = ecoli_core_reaction_kcats[rid] * 3.6, #src + kcat_forward = kcat, #src + kcat_reverse = kcat, #src ) #src end end gene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses) -Symbolics.@variables capacitylimitation -parameter_values[capacitylimitation] = 50.0 # mg enzyme/gDW +@variables capacitylimitation +parameter_values[:capacitylimitation] = 50.0 # mg enzyme/gDW km = build_kinetic_model( model; @@ -100,15 +102,15 @@ ec_solution_cobrexa = enzyme_constrained_flux_balance_analysis( #src @test isapprox(ec_solution.objective, ec_solution_cobrexa.objective; atol = TEST_TOLERANCE) #src # This solution contains many inactive reactions -sort(abs.(collect(values(ec_solution.fluxes)))) +sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last)) -@test any(abs.(collect(values(ec_solution.fluxes))) .≈ 0) #src +@test any(values(ec_solution.fluxes) .≈ 0) #src # And also many inactive gene products. -sort(abs.(collect(values(ec_solution.gene_product_amounts)))) +sort(collect(ec_solution.gene_product_amounts), by=last) -@test any(round.(abs.(collect(values(ec_solution.gene_product_amounts))); digits = 8) .≈ 0) #src +@test any(isapprox.(values(ec_solution.gene_product_amounts), 0, atol=1e-8)) #src # With theory, you can show that this introduces flux variability into the # solution, making it non-unique, and consequently non-differentiable. To fix @@ -121,7 +123,7 @@ flux_zero_tol = 1e-6 # these bounds make a real difference! gene_zero_tol = 1e-6 pruned_reaction_isozymes = - prune_reaction_isozymes(reaction_isozymes, ec_solution, gene_zero_tol) + prune_reaction_isozymes(reaction_isozymes, ec_solution, flux_zero_tol) pruned_gene_product_molar_masses = prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol) @@ -147,7 +149,7 @@ ec_solution2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_ ec_solution2 # no zero fluxes -sort(abs.(collect(values(ec_solution2.fluxes)))) +sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last)) # no zero genes sort(abs.(collect(values(ec_solution2.gene_product_amounts)))) @@ -159,25 +161,21 @@ sort(abs.(collect(values(ec_solution2.gene_product_amounts)))) k in intersect(keys(ec_solution.fluxes), keys(ec_solution2.fluxes)) #src ) #src -# prune the kcats, leaving only those that are actually used -pruned_kcats = [ - begin - x = first(values(v)) - isnothing(x.kcat_forward) ? x.kcat_reverse : x.kcat_forward - end for v in values(pruned_reaction_isozymes) -] - -parameters = [capacitylimitation; kcats] +parameters = [:capacitylimitation; kcats] -sens, vids = differentiate( +pkm_kkt, vids = differentiate_prepare_kkt( pkm, pkm.objective.value, + parameters +) + +sens = differentiate_solution( + pkm_kkt, x_vals, eq_dual_vals, ineq_dual_vals, parameter_values, - parameters; scale = true, # unitless sensitivities ) diff --git a/dev/2-differentiate-enzyme-model/index.html b/dev/2-differentiate-enzyme-model/index.html index 4224b4b..cd86b99 100644 --- a/dev/2-differentiate-enzyme-model/index.html +++ b/dev/2-differentiate-enzyme-model/index.html @@ -1,13 +1,13 @@ Differentiating enzyme constrained metabolic models · DifferentiableMetabolism.jl
GitHub

Differentiating enzyme constrained metabolic models

using DifferentiableMetabolism
 using AbstractFBCModels
-using Symbolics
 using ConstraintTrees
 using COBREXA
 using Tulip
 using JSONFBCModels
 import Downloads: download
 using CairoMakie
+using FastDifferentiation
 
 !isfile("e_coli_core.json") &&
     download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
@@ -19,12 +19,10 @@
 model.reactions[glc_idx]["lower_bound"] = -1000.0
-1000.0

constrain PFL to zero

pfl_idx = first(indexin(["PFL"], rids))
 model.reactions[pfl_idx]["upper_bound"] = 0.0
 
-kcats = vcat([Symbolics.@variables $x for x in Symbol.(keys(ecoli_core_reaction_kcats))]...)
+rid_kcat = Dict(k => FastDifferentiation.Node(Symbol(k)) for (k,_) in ecoli_core_reaction_kcats)
+kcats = Symbol.(keys(ecoli_core_reaction_kcats))
 
-rid_kcat = Dict(zip(keys(ecoli_core_reaction_kcats), kcats))
-
-parameter_values =
-    Dict(kid => ecoli_core_reaction_kcats[rid] * 3.6 for (rid, kid) in rid_kcat) # change unit to k/h
+parameter_values = Dict{Symbol, Float64}()
 
 reaction_isozymes = Dict{String,Dict{String,ParameterIsozyme}}() # a mapping from reaction IDs to isozyme IDs to isozyme structs.
 for rid in AbstractFBCModels.reactions(model)
@@ -32,8 +30,12 @@
     isnothing(grrs) && continue # skip if no grr available
     haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available
     for (i, grr) in enumerate(grrs)
+
+        kcat = ecoli_core_reaction_kcats[rid] * 3.6 # change unit to k/h
+        parameter_values[Symbol(rid)] = kcat
+
         d = get!(reaction_isozymes, rid, Dict{String,ParameterIsozyme}())
-        d["isozyme_"*string(i)] = ParameterIsozyme(
+        d["isozyme_$i"] = ParameterIsozyme(
             gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes
             kcat_forward = rid_kcat[rid],
             kcat_reverse = rid_kcat[rid],
@@ -43,8 +45,8 @@
 
 gene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses)
 
-Symbolics.@variables capacitylimitation
-parameter_values[capacitylimitation] = 50.0 # mg enzyme/gDW
+@variables capacitylimitation
+parameter_values[:capacitylimitation] = 50.0 # mg enzyme/gDW
 
 km = build_kinetic_model(
     model;
@@ -74,51 +76,51 @@
   :isozyme_flux_reverse_balance => ConstraintTrees.Tree{Float64}(#= 32 elements =#)
   :isozyme_forward_amounts      => ConstraintTrees.Tree{Float64}(#= 60 elements =#)
   :isozyme_reverse_amounts      => ConstraintTrees.Tree{Float64}(#= 32 elements =#)
-  :objective                    => 0.750518

This solution contains many inactive reactions

sort(abs.(collect(values(ec_solution.fluxes))))
95-element Vector{Float64}:
-   0.0
-   0.0
-   0.0
-   0.0
-   0.0
-   0.0
-   0.0
-   0.0
-   0.0
-   4.685845906545776e-15
-   ⋮
-  33.89684016777724
-  33.89684016777725
-  38.668072869856005
-  38.668072869856005
-  39.54058803115938
-  64.6666449207555
-  64.66664492281971
- 101.83896464908304
- 101.83896465010721

And also many inactive gene products.

sort(abs.(collect(values(ec_solution.gene_product_amounts))))
124-element Vector{Float64}:
- 8.141678306020198e-17
- 8.141678306020198e-17
- 8.141678306020198e-17
- 9.846484300151926e-13
- 1.0446152588779501e-12
- 1.0446152588779501e-12
- 1.097789630627885e-12
- 1.256902936168616e-12
- 1.295043723786085e-12
- 2.807358103195605e-12
- ⋮
- 0.022586033326179002
- 0.022586033355673742
- 0.022586033355673742
- 0.04348653861416314
- 0.0731266615462882
- 0.11726698379336735
- 0.11726698379336735
- 0.12079853655910855
- 0.16335511688953944

With theory, you can show that this introduces flux variability into the solution, making it non-unique, and consequently non-differentiable. To fix this, one need to prune the model, to include only the active reactions and genes. This can be differentiated uniquely. Below we build a pruned kinetic model, by removing all the reactions, metabolites, and genes that are not active.

flux_zero_tol = 1e-6 # these bounds make a real difference!
+  :objective                    => 0.750518

This solution contains many inactive reactions

sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))
95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
+    :EX_fru_e => 0.0
+    :EX_fum_e => 0.0
+ :EX_gln__L_e => 0.0
+ :EX_mal__L_e => 0.0
+     :FRUpts2 => 0.0
+     :FUMt2_2 => 0.0
+      :GLNabc => 0.0
+     :MALt2_2 => 0.0
+         :PFL => 0.0
+    :EX_for_e => 4.685845906545776e-15
+              ⋮
+         :PGK => -33.89684016777724
+        :GAPD => 33.89684016777725
+    :EX_h2o_e => 38.668072869856005
+        :H2Ot => -38.668072869856005
+      :EX_h_e => 39.54058803115938
+      :NADH16 => 64.6666449207555
+       :CYTBD => 64.66664492281971
+    :SUCCt2_2 => 101.83896464908304
+      :SUCCt3 => 101.83896465010721

And also many inactive gene products.

sort(collect(ec_solution.gene_product_amounts), by=last)
124-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
+ :b0809 => 8.141678306020198e-17
+ :b0810 => 8.141678306020198e-17
+ :b0811 => 8.141678306020198e-17
+ :b4014 => 9.846484300151926e-13
+ :b0726 => 1.0446152588779501e-12
+ :b0727 => 1.0446152588779501e-12
+ :b1702 => 1.097789630627885e-12
+ :b1479 => 1.256902936168616e-12
+ :b2976 => 1.295043723786085e-12
+ :b3403 => 2.807358103195605e-12
+        ⋮
+ :b2417 => 0.022586033326179002
+ :b2415 => 0.022586033355673742
+ :b2416 => 0.022586033355673742
+ :b2779 => 0.04348653861416314
+ :b1779 => 0.0731266615462882
+ :b0978 => 0.11726698379336735
+ :b0979 => 0.11726698379336735
+ :b3528 => 0.12079853655910855
+ :b2926 => 0.16335511688953944

With theory, you can show that this introduces flux variability into the solution, making it non-unique, and consequently non-differentiable. To fix this, one need to prune the model, to include only the active reactions and genes. This can be differentiated uniquely. Below we build a pruned kinetic model, by removing all the reactions, metabolites, and genes that are not active.

flux_zero_tol = 1e-6 # these bounds make a real difference!
 gene_zero_tol = 1e-6
 
 pruned_reaction_isozymes =
-    prune_reaction_isozymes(reaction_isozymes, ec_solution, gene_zero_tol)
+    prune_reaction_isozymes(reaction_isozymes, ec_solution, flux_zero_tol)
 
 pruned_gene_product_molar_masses =
     prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol)
@@ -149,65 +151,67 @@
   :isozyme_flux_reverse_balance => ConstraintTrees.Tree{Float64}(#= 5 elements =#)
   :isozyme_forward_amounts      => ConstraintTrees.Tree{Float64}(#= 26 elements =#)
   :isozyme_reverse_amounts      => ConstraintTrees.Tree{Float64}(#= 5 elements =#)
-  :objective                    => 0.750518

no zero fluxes

sort(abs.(collect(values(ec_solution2.fluxes))))
48-element Vector{Float64}:
-   0.19190755046855745
-   0.7505183827479407
-   0.8097342831462997
-   0.8097342831463148
-   0.8097342831463405
-   0.8097342831463762
-   1.739538998893496
-   2.010476135065505
-   2.010476135065506
-   2.150685477602044
-   ⋮
-  33.89684011077103
-  33.89684011077104
-  38.66807309949494
-  38.66807309949494
-  39.54058824535901
-  64.66664537982322
-  64.66664537982322
- 101.83896538885946
- 101.83896538885948

no zero genes

sort(abs.(collect(values(ec_solution2.gene_product_amounts))))
47-element Vector{Float64}:
- 0.0005876895265927602
- 0.000593428174368933
- 0.0008428927106180381
- 0.0019854019751287073
- 0.0025855775424920503
- 0.002828557018354045
- 0.0030325373898386107
- 0.003279421656136242
- 0.003609307833900057
- 0.004477678876501782
- ⋮
- 0.022586033324645538
- 0.022586033324645538
- 0.022586033324645538
- 0.04348653850566591
- 0.07312666138287223
- 0.11726698687786197
- 0.11726698687786197
- 0.12079853743661037
- 0.16335511657978186

prune the kcats, leaving only those that are actually used

pruned_kcats = [
-    begin
-        x = first(values(v))
-        isnothing(x.kcat_forward) ? x.kcat_reverse : x.kcat_forward
-    end for v in values(pruned_reaction_isozymes)
-]
+  :objective                    => 0.750518

no zero fluxes

sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))
95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
+    :EX_fru_e => 0.0
+    :EX_fum_e => 0.0
+ :EX_gln__L_e => 0.0
+ :EX_mal__L_e => 0.0
+     :FRUpts2 => 0.0
+     :FUMt2_2 => 0.0
+      :GLNabc => 0.0
+     :MALt2_2 => 0.0
+         :PFL => 0.0
+    :EX_for_e => 4.685845906545776e-15
+              ⋮
+         :PGK => -33.89684016777724
+        :GAPD => 33.89684016777725
+    :EX_h2o_e => 38.668072869856005
+        :H2Ot => -38.668072869856005
+      :EX_h_e => 39.54058803115938
+      :NADH16 => 64.6666449207555
+       :CYTBD => 64.66664492281971
+    :SUCCt2_2 => 101.83896464908304
+      :SUCCt3 => 101.83896465010721

no zero genes

sort(abs.(collect(values(ec_solution2.gene_product_amounts))))
+
 
-parameters = [capacitylimitation; kcats]
 
-sens, vids = differentiate(
+
+parameters = [:capacitylimitation; kcats]
+
+pkm_kkt, vids = differentiate_prepare_kkt(
     pkm,
     pkm.objective.value,
+    parameters
+)
+
+sens = differentiate_solution(
+    pkm_kkt,
     x_vals,
     eq_dual_vals,
     ineq_dual_vals,
     parameter_values,
-    parameters;
     scale = true, # unitless sensitivities
-)
([0.9447742842823873 -0.0 … -0.0 -0.0; 1.1401065257754806 0.0 … 0.0 0.0; … ; 1.140106525774991 0.0 … 0.0 0.0; 0.9447742842823852 0.0 … 0.0 0.0], Tuple{Symbol, Symbol, Vararg{Symbol}}[(:fluxes, :ACKr), (:fluxes, :ACONTa), (:fluxes, :ACONTb), (:fluxes, :ACt2r), (:fluxes, :ATPM), (:fluxes, :BIOMASS_Ecoli_core_w_GAM), (:fluxes, :CO2t), (:fluxes, :CS), (:fluxes, :CYTBD), (:fluxes, :ENO)  …  (:isozyme_forward_amounts, :ENO, :isozyme_1), (:isozyme_forward_amounts, :CYTBD, :isozyme_2), (:isozyme_forward_amounts, :CS, :isozyme_1), (:isozyme_forward_amounts, :ACONTb, :isozyme_1), (:isozyme_forward_amounts, :ACONTa, :isozyme_1), (:isozyme_reverse_amounts, :RPI, :isozyme_2), (:isozyme_reverse_amounts, :PGM, :isozyme_3), (:isozyme_reverse_amounts, :PGK, :isozyme_1), (:isozyme_reverse_amounts, :GLUDy, :isozyme_1), (:isozyme_reverse_amounts, :ACKr, :isozyme_1)])

look at glycolysis and oxidative phosphorylation only

subset_ids = [
+)
127×63 Matrix{Float64}:
+ 0.944774     -0.0  0.00395506    0.000578273  -0.0  -0.0   0.126918     …  -0.0   0.00340805   0.0375147    -0.0  -0.0  -0.0  -0.0  -0.0
+ 1.14011       0.0  0.00477276    0.000697831   0.0   0.0   0.153158         0.0   0.00411266   0.0452708     0.0   0.0   0.0   0.0   0.0
+ 1.14011      -0.0  0.00477276    0.000697831  -0.0  -0.0   0.153158        -0.0   0.00411266   0.0452708    -0.0  -0.0  -0.0  -0.0  -0.0
+ 0.944774     -0.0  0.00395506    0.000578273  -0.0  -0.0   0.126918        -0.0   0.00340805   0.0375147    -0.0  -0.0  -0.0  -0.0  -0.0
+ 4.76087e-13  -0.0  2.76756e-11  -1.02542e-11  -0.0  -0.0   1.59896e-11     -0.0  -6.04713e-11  2.76927e-11  -0.0  -0.0  -0.0  -0.0  -0.0
+ 1.14011      -0.0  0.00477276    0.000697831  -0.0  -0.0   0.153158     …  -0.0   0.00411266   0.0452708    -0.0  -0.0  -0.0  -0.0  -0.0
+ 0.996052      0.0  0.00416972    0.000609659   0.0   0.0   0.133806         0.0   0.00359302   0.0395508     0.0   0.0   0.0   0.0   0.0
+ 1.14011       0.0  0.00477276    0.000697831   0.0   0.0   0.153158         0.0   0.00411266   0.0452708     0.0   0.0   0.0   0.0   0.0
+ 0.992186      0.0  0.00415353    0.000607293   0.0   0.0   0.133287         0.0   0.00357908   0.0393973     0.0   0.0   0.0   0.0   0.0
+ 0.994176     -0.0  0.00416186    0.00060851   -0.0  -0.0   0.133554        -0.0   0.00358625   0.0394763    -0.0  -0.0  -0.0  -0.0  -0.0
+ ⋮                                                    ⋮                  ⋱   ⋮                                            ⋮          
+ 0.992186      0.0  0.00415353    0.000607293   0.0   0.0   0.133287         0.0   0.00357908   0.0393973     0.0   0.0   0.0   0.0   0.0
+ 1.14011      -0.0  0.00477276    0.000697831  -0.0  -0.0   0.153158        -0.0   0.00411266   0.0452708    -0.0  -0.0  -0.0  -0.0  -0.0
+ 1.14011      -0.0  0.00477276    0.000697831  -0.0  -0.0   0.153158     …  -0.0   0.00411266   0.0452708    -0.0  -0.0  -0.0  -0.0  -0.0
+ 1.14011       0.0  0.00477276    0.000697831   0.0   0.0   0.153158         0.0   0.00411266   0.0452708     0.0   0.0   0.0   0.0   0.0
+ 1.14011      -0.0  0.00477276    0.000697831  -0.0  -0.0   0.153158        -0.0   0.00411266   0.0452708    -0.0  -0.0  -0.0  -0.0  -0.0
+ 0.994176     -0.0  0.00416186    0.00060851   -0.0  -0.0   0.133554        -0.0   0.00358625   0.0394763    -0.0  -0.0  -0.0  -0.0  -0.0
+ 0.999009     -0.0  0.0041821     0.000611469  -0.0  -0.0  -0.865796        -0.0   0.00360369   0.0396682    -0.0  -0.0  -0.0  -0.0  -0.0
+ 1.14011       0.0  0.00477276    0.000697831   0.0   0.0   0.153158     …   0.0   0.00411266   0.0452708     0.0   0.0   0.0   0.0   0.0
+ 0.944774      0.0  0.00395506    0.000578273   0.0   0.0   0.126918         0.0   0.00340805   0.0375147     0.0   0.0   0.0   0.0   0.0

look at glycolysis and oxidative phosphorylation only

subset_ids = [
     r["id"] for r in model.reactions[indexin(string.(keys(pkm.fluxes)), rids)] if
     r["subsystem"] in ["Glycolysis/Gluconeogenesis", "Oxidative Phosphorylation"]
 ]
@@ -219,24 +223,18 @@
 iso_ids = [v[2] for v in vids[iso_idxs]]
 
 param_idxs = findall(x -> string(x) in subset_ids, parameters)
-param_ids = parameters[param_idxs]

\[ \begin{equation} -\left[ -\begin{array}{c} -\mathtt{PGK} \\ -\mathtt{PDH} \\ -\mathtt{PGM} \\ -\mathtt{CYTBD} \\ -\mathtt{FBA} \\ -\mathtt{GAPD} \\ -\mathtt{NADH16} \\ -\mathtt{PFK} \\ -\mathtt{PGI} \\ -\mathtt{TPI} \\ -\mathtt{ENO} \\ -\end{array} -\right] -\end{equation} - \]

Flux sensitivities

f, a, hm = heatmap(
+param_ids = parameters[param_idxs]
11-element Vector{Symbol}:
+ :PGK
+ :PDH
+ :PGM
+ :CYTBD
+ :FBA
+ :GAPD
+ :NADH16
+ :PFK
+ :PGI
+ :TPI
+ :ENO

Flux sensitivities

f, a, hm = heatmap(
     sens[flux_idxs, param_idxs]';
     axis = (
         xlabel = "kcat",
@@ -260,4 +258,4 @@
     ),
 )
 Colorbar(f[1, 2], hm)
-f
Example block output

This page was generated using Literate.jl.

+fExample block output

This page was generated using Literate.jl.

diff --git a/dev/3-parameter-estimation.jl b/dev/3-parameter-estimation.jl index ad3136c..63727e0 100644 --- a/dev/3-parameter-estimation.jl +++ b/dev/3-parameter-estimation.jl @@ -16,7 +16,7 @@ using DifferentiableMetabolism using AbstractFBCModels -using Symbolics +using FastDifferentiation using ConstraintTrees using COBREXA using Clarabel @@ -39,7 +39,7 @@ model.reactions["r2"].lower_bound = -1000.0 # ![simple_model](./assets/simple_model_pruned.svg) -kcats = Symbolics.@variables r3 r4 +@variables r3 r4 reaction_isozymes = Dict( "r3" => Dict( @@ -60,9 +60,9 @@ reaction_isozymes = Dict( gene_product_molar_masses = Dict("g1" => 20.0, "g2" => 10.0) -Symbolics.@variables capacitylimitation +@variables capacitylimitation -true_parameter_values = Dict(capacitylimitation => 50.0, r3 => 2.0, r4 => 3.0) +true_parameter_values = Dict(:capacitylimitation => 50.0, :r3 => 2.0, :r4 => 3.0) km = build_kinetic_model( model; @@ -103,12 +103,17 @@ km *= bound = nothing, ) -estimated_parameters = Dict(capacitylimitation => 50.0, r3 => 5.0, r4 => 1.0) # initial values -parameters = [r3, r4] # will differentiate against these two parameters +estimated_parameters = Dict(:capacitylimitation => 50.0, :r3 => 5.0, :r4 => 1.0) # initial values η = 0.1 # learning rate losses = Float64[] +kmKKT, vids = differentiate_prepare_kkt( + km, + km.loss.value, + [:r3, :r4, :capacitylimitation], +) + for k = 1:150 sol2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters( @@ -122,16 +127,14 @@ for k = 1:150 push!(losses, sol2.loss) - sens, vids = differentiate( - km, - km.loss.value, + sens = differentiate_solution( + kmKKT, x_vals, eq_dual_vals, ineq_dual_vals, estimated_parameters, - parameters; ) - measured_idxs = [1, 3, 11, 12] + measured_idxs = [1, 3, 12, 11] x = [ sol2.fluxes.r1, @@ -143,8 +146,8 @@ for k = 1:150 dL_dx = x - measured # derivative of loss function with respect to optimization variables dL_dkcats = sens[measured_idxs, :]' * dL_dx # derivative of loss function with respect to parameters - estimated_parameters[r3] -= η * dL_dkcats[1] - estimated_parameters[r4] -= η * dL_dkcats[2] + estimated_parameters[:r3] -= η * dL_dkcats[1] + estimated_parameters[:r4] -= η * dL_dkcats[2] end lines(losses; axis = (xlabel = "Iterations", ylabel = "L2 loss")) @@ -156,7 +159,6 @@ estimated_parameters true_parameter_values -@test abs(estimated_parameters[r3] - true_parameter_values[r3]) <= 0.1 #src -@test abs(estimated_parameters[r4] - true_parameter_values[r4]) <= 0.1 #src -@test abs(estimated_parameters[r4] - true_parameter_values[r4]) <= 0.1 #src +@test abs(estimated_parameters[:r3] - true_parameter_values[:r3]) <= 0.1 #src +@test abs(estimated_parameters[:r4] - true_parameter_values[:r4]) <= 0.1 #src @test all(losses[2:end] .<= losses[1:end-1]) #src diff --git a/dev/3-parameter-estimation/98c59523.png b/dev/3-parameter-estimation/98c59523.png new file mode 100644 index 0000000000000000000000000000000000000000..4f80929e8819bbc6436445df69b2a79eb23fc5a9 GIT binary patch literal 32477 zcmd43cT`l_w=RmIr38Uel%OD#B2j|m3=N8sRY0(P5Cy=$+v=bGW0-<)fe&u^++CnGsaLPSJFCV%7V zZ6czhIYdN9@R-AJN0P!I6aG4OPx1OyqJ8AQx0NZOL_~~4@>l;>e-t}6;`%sj!c}N* zyII7n-}I?b=rwgJ>eJT*pFDqIn#6>Yvzp=VO|uiQC}`;|h~Be|?oG~AwmLi3`{Cu4 z=N*wH-H}?yzkE3UNne`m-Pd=#U#rNDe9oBK{r)lPkZGvN`ANz6HiIegsVEWZfbY!_ zvt^57VNDZ=V6`bs1my1T25pI=j8JW52w zdiyYrh)6o*6y^{S(TmsrzZk+E>Uc=vE_;|l83@RzQ>wS{)AJ!cx{M^*N9J_B-J@}X zgM(z}1Xl=MW?zGGBl`;<6K-p0^b8M+Km76a)v5g(#1cKdnY!ue=^XIG&i*RvFJbA` z)m5v$EKW{NdO>Tg&ceA~vaV;fr!dsqcPq?0-t`%HSLWva`a-E1EBxTwm*O&IJ=1bM9_dS=CY?`@84^%mC!(1za708Kvk7sm<;2zdT{!r-wv5h06V ziT7B73HYL#m)oTF2d%Otc7A@&u2^bR%_iA6vW>-rR&VqGbopfc~8C|cUmJwd8SgzWc;HJJ~ zoP0YK1}tRcx02w**ZS|v_YxBk1;<7nF;Y|O>P%9U#l&xXVOjHDoy%_{$5rfY&4Lz2 z^l6ERKAvQp6u{5+W}XwW)x(Uhbg7Q*_tgA-lHgdT+XT)c%=?87b3K427jz?Nkd~bt zFXm#4F&ZqeKs@GV>^*Y)wz;+SaFyqFvT{tE=i0!KZMjXyyKA0XE1Q!-)lrW6u*wG3 zchVbQF9zw^fBMTeMJ+{Tx7DS}p^nzfam(j^ec=2`oXe=|)Y}(gjOEVrFt^&`Xyv)H zXUq0Sd}K{V%A8tSS}q|!SF&fHmS|{d+Kp9J78MmGCo>at9v|s$Z_n}C-Le)O8Xk^y z^kEp2jg zGUVcQYF zZEWVZR%gNE;~w2RTnkRh#5qYm-kqwpOtDbDShZ{5Ii29QwB4W*E1dhky-HP-@ZO>; z*%zBty|)n@|579z!6m%=%9 zTP3&WC_Mf>0wMO0Jxr(K(Jx9>@sOs8yn(q9m$4Co2m)=lUWj0yVpU+ufqXOXz2$g= zJmZFJ@7>v*S2Y#+t)jD1`^!=-{@~=-IF_VoWH$aKmfGl4`=`ICne;(zA44upDx1K|rXT-*22Eah>W>L4{heg{ zQ!M_@-9jzmn;+cw9zE=wfO)U;`|z+kDYs}0Is0OZdIXvENr{N;-}CXF4O41Z;5NVS z;E)%FJSPdyMZHqr`&(9xYd`_O{$ji|O_;Tt3!}9)6~Ovq6%76sn%F<$gdh80`DbXeYs5Jy`ng*+@&xQ!Yg$Q zE`n4!8tfb^2-23_6yiLIX_62=!-I2aown_a5tHUH%%52EBexqTfbjr+=(iN>m=rWDCfbn z&!uT%Wc+e+<&iWi3L}|y`%To*IfLONC8LO2s09!~8V(e85gxK4Lv$o?*vOifvpnD| zEw~yoxb7=>Vf!_>9e~Qm^j~k0(AOM8^5a`s_`EY!;$mViVgP#3Mpxd`o2UTDfV0Rm zh-Qp4kmZ!gfkKMDpfg|kZ!_m2%0|^eg9v~=$XvN?HP~*9dcreMSxh51pFT{RY(IB3 zQQ(i=PWNI3`%MscK~^~zCx914KB|Ye9^6JOr(KRN0@nM#mCPB_sD&yPF;i1AuQR&z z|F+KmVd?+5#T+AuS`*ZM4%G#^P6cOg<#W;KY7s5|zZB z&QoG8s`r~0KD;Ln9ZHSCpKLtsOk4OA z{NEg~Bq7`ZGml7Bpvo9<)4A=e465KGFk%iB))#puL6Zkzfr#}#-1-9_h#RoRjP9DCPKLP1 z*JNIZ4P+6#e-w#Bjr-qjB`B5g@Us7LfY+RW%m>mVDrG;@c_|eu!$lpf^b&CqEarQ6 zUdIPi<+6+cYE5#9DsiCFK^(TOK@Nerl%qAmg`zQc5m*>x#@40w46LA17ahkxIFLhu zG2p-}AYH`fS)k^ogkO0tl~?DHp(aA+lShs5TAGFh1p~+w-wzal+(OrZybE>sSL5BuN)1);2fe3ea)?|I|x9ePohZ1D_`RH!NP3JV z@yTIh`SaG-^0UaV$P=&8fS3lESx|+bx^V}%6a6CIsMgTI2zV zdKG?t)X>P(|M80yM|t)z5*$8xQcr~makmRtIdDjjy8dX5xd@t{Bc*9L&APxg^uQ{@ za+aXvgcz5@6I$7O;ss8m>IUdN^XZ=)@AxIL1NAIAXjuH&DIRlO4{b zUuN18{XIx(uUk!ZcO^y2bLUq*i+-QcfsyY+Nhe^{>#^VvzrES!eYYgf@Y6rrz*pn< z2=P!%m1TR*TSY0wdLKu_BCs~fSEu~j3y#UwH|J!8r(zvoasiKxpPB}qtJ`DydxM2m zb|vFa^|)LXBhlcd%UpPxR&=f9?W;`ZKb5*l$-V=pYXa2|2$gW*=A6|Va` zZTlX3n^US$yq7w7DM-Mo*x1<8?0!HI>n>)V*RLWq1x12``3rX;sGSO3 zmnZL5x)qA74)))S@XoRQob1B>V!GO5`T=wSBAO-?6cmoiOqZsmr(0QC!T=35wNOXS zkB2PH&HHnVWCdevY`NYk)B1B@1^M}h-DWfQe^Tt{z^vFBlL&q1E@dI7ty}l&0to;1 zc5s;LWWU%v%I)Y1ou$I5yp>U2Ds0Cmk z=PwtGklm;_)bY@h@WXPI{QFOgZJEW&bcZ~34}BgmcX~14(*>87q)Qy8T7x|2?)y`7 zFZ5-Py>S`cyzaenU0$=>Q^yJ8Z+`1CM6*)}aZQX})h5)h4W{`p!gFkUG{RA2ggxKv z$}5!RBuDlkEtC<^JPvP_Rhr!2X=pCNYl)KK@JKN{^ig)wG}2`-;9pzW`&-#-wqnq$ zky;7~a=foU^4`R_-~hUwMjjc@Gih3)u&v%*BeyNz_~xH@0za@47jg?lZO$t}0nQ#? zhASlWGDtVfT&>yXCseg=S)^V|q47-mgg@ci?342ql!CJTf@_Qccj{2Q^#kNGIVfEK znbO;F1t?Dnm8Qwe_``zD3o~%mNK5UpAO_!IDdPNY01ee+fMEYcl*2{^pe3{}BZJ=| z@pPcC(1ln0fpnXbn8Y65-H`hjq4-{2*m~H;h{mcy1PAD+{)}%6s%zKbHKoK0}Q~k$L zg(8}O!Y2VOd=Q*YhfIT(9>^kZrv`J$`(huKwi{bQit7)Orr|q)e=O$yXK_Wf}Aa|pu*;#Q}G5*FkgphT1b ztw9#fn4ps(Mv~VyRDR(j`A=^mkM#p-&cAcuU;PpAy9}ds0m5)|7<_n4!BB{*KYD=- ze~9QD&>d}H%Kup#FF8_N41ze#sOKDlfoL`&$KbDYAw>%C$GQwGCgvehU*~Yc=TQgs zHIyc%0|yx!c=iHgmKV(C=ePA&2IM~rwLo56xb@BQsQJ%<(`9i!FW`_S?8WAGCG zblaVE3UMa-8@yO9tsF+YeUFBwzGKHV=<(r$n`#F>ZME zgMB{rz`J2w)#hXj4&M;ULY-ByHIuGeWPKc0vA?^>ttSp$I})Mgtb;--*qHG$2L~yx zyRVNrOMGwRI}YF2*m#_h`83xgMLg?e-uJV=(ZcyKYR;JN5Nrdz^$@5?h_XXyn6 z{CFfcdz=-Y`H)z(Ub&dmMfgH8pPP43nM4@-@!)xYpusbh=qtJCpnVFP1+W*e%pDdM zrvGT%S1nnobg7ZUkDNZ*`i+Eoa!%6+v@x{HM2C^Z5FFjW7ebkp2dB1DQ{uM%eL}%} z`!>IhGc8)YtvURNn`rOuel2DH&yjF8`J~%8e0z6n({m)YMqg3jfjvRLd?7JNQl})( zZFNRF-*k1hFT0|`1Df3VW^F1t6>(hm6OsNA{c+-&&E=^eaB0|5fOfdqo&q$NeW~vv zjq~oF9sz!S@%3SceAAY)^-=fG%>1c`f)bEd0Ffbe$47GvBgS|NLK3KLd78>%e^pYgBkltG;(tE4CM2i{oSZ=}A{AdUz z)kL;CICYEg*wOvHb!x8jF2fIeV_b{T0GOBO#fFE4`C@_4kZPXy-uZ=>n0h(~rDAHDBPfWt z|6h7^G@AV4l>S-tg#ND*Y!JadM}&cw>-e{l2j5$Pn_z2Tv&F167TIjT*7WDzm)dMl zWxivs+B-vswo44T<(+{mcY{Z4!6r_zy~&qnCy`*Xu&}@wEl;*=@;~anbsh0^#gwk+ zX9hzLTCvVd1Gp1NE#Lta19$b#608Q<4geSe{q885n3z!CeKh`>O)U#Hm=a|D0P)(w zkW7Ts_iG{cj3;pKy9kIYNAp|6&qBL><_U#l^>VI6Uk~WNx$j z-w@Q-I^cf&im!Zrs_y+zg&OCs>qHnwR>`~MWYhg~g4X@^E%gxCCA#83MAoRESjo5lxdX^~ILhj|+p7!>iP>M zH)(!jZiu}fkKNLbhE0bl*w~HWHi&C6tq-8>g5u}>>QA=FE%6R$1hYykb7PmV! zuGXUDDGq;KPGlJ9pt1;30~m6T^x-nn(g3*9$u~Y%Z?@Q}MDmPdd)1hD2{H#AlC@4x zpBzwQYpe+I&bYzuj3$NAmD-EC_t1)@Xn__?za!^whPL*CPI$ z1;F@`opY|IcX1pV3IUD`h%X@+_FS(im-4{c(3$LeICmB8K^k%&Hu~|0J5?UXhgi`b z^kMs<6;%4C$LbD4B9w=aE<4g9D*Zz>p$aJd57iX%;dzKM#x>(pu6@upK-3CJMGNxC zah6yfsAm6ua;(sG3U<}Jwz}eYNatk$36cQsmWg4xfTutP6%Y{EOi5K|k7(-Y z?%o+0+Z}<;HQ2sxeRWDL!|!%_$iXH^rK8K%*3Oia0o3I!rid>A@IdO5*q+O0*!+53 z1|zjIPEox(Z|CA5(ZPo{qf>N`i~|@#Z^hs#wrYE1(QDqy;GBe;qhO3@oBJrA``Vl$ zS1A+T9>nF=Ey@s6d~sJ0EojVji;vX+e2d`iVHecEYjeW=%&prGx95w4&P$XP6|JQk z>{_~npRR0|f_rjZQBV6##4p$K83{p%oVI^;P@6=59yCy-Jhe*Oa2&n-e@ zu_`H3ZSmeiHpNb6MkkvN`fWC|z0$(2OFz`v3m;B2R&T~sTjUc9TXep!-te!+;5;Tc z3|wbYQ_Av-v_#Nuq$BTo4`kX+kCd<950Vgc zn&}31!R<*)k~*4vy^wmug)N4W}l~;+v)e>xd1oNHha_aiR{TSva-o&tFp{E*1|z`bo(qP;;1No zQ9)Mr^EGOOgHv;YxV?X+p2h@i${!CC8qm$1#*p(!oq|}Fn8Cvy@Ls;JQ_ZC7 zfI8W~Czc_AtHZrXa8}0n&RItI+hT;vg504kB1;OL^Urh%)iZ;-i?fYHHavrVO zSZ!U{V_*DWr&M%`i0I#EFm^GrQF4>m&L&m);MOZbk&8~T?MPqs-m1&?EvQTU5+mR7 zT1t=YJP{IF{5=qG%DyR_vv|zQy||M8uf$^u3)reuG@0Ks=J!fGhHm28OjnVQ164VD z_ADS?7hwf=Ir#&0@iY<5Npp-)1LciUDz2h;g(QS-*uDtq*dXIr94Y3fm+7y_HyX-1 zV%G_kOvp*h_Xq{w`{>Rql&zOZd_7I>*zyX%{?HsVmdf`Iiz)^I0o%hyx6GXcE76ix z$R~$Y_{BkKtaVzB`tU)eO98iG!!T`Ym<1Or`J8Zp3*zXM?lknNr~|?7pGcc_{qaW(j72y6D!*FaL3?llX3#u z2UcLPw3Z%ZE;}1EE^Ad_gzrHE2E^k^Bu1Mxh~gtU4m1JlS9Hy8{#Z4Dt;wpJbP_;r zKYj;Kya|6PClm$SlfVBj5mSs<>Sp)Rns@gNp!LV7%6M6{o~+;V^mt4R{WaJiy9Oo^ zVvzWM`{+2n_5%-ZnWL;SdaU6% zGVp!k_7p?2_wSPq=&jcwyW3&OvKS7*`HaCL4#Gp2*%Udx&nEHmooHm$a0renm2jVUf}vR#2UR(8>dn79s-B!xIRcOMsxETt3fy|% z1g*_ejPPVt5~1PZwjUHVKf7aeUvWV{=7w~~Uo`-qS^DKgaB_$s^QzC1Ho+72Q{U2@ z(P!Klo{|)Tx)Ye|CH)&qt+Y0@(3?43_>BY)c~EM5pzT%Gw;K=EB7H8Us$<4$(`>|M zV$g}FZLgiioI7_8IV3heRLr9PHRPUj{M?2Su*9dylV+oN0pMjXRG0~e$;*pJ-K?Pg zKzI{{*Uu-=Yk1!zlBdML_cKhO2Zwj%5cK?s|Acy&WBSICai>20^{)6megb>o&8Yxd z$EF>jCrB{ed<}@0jMZvrDzJY>ig}H#q$qrh!Zk zg_`CRgVg4a(8ESnpE=xu=)3JUlj+4~W2n?^#Fk4A$!kKaI5N;D-0ORIHX#+tTUK6Z zCsT}8@6A>_yQu=16;jU&g427DH|px@fI-mzWQh{3Xrqw9;KYB`l}W|${=kMEs~`o9 z!Qn6HSzo%e7eelq25NyrSWwEmCp|>2=jrJQfG8LPckMb09T|Lb91xNwR%r$&+D2|R zNFq&seSLjR&4>oNjlN&qhswPVdCs(C#uQsOz?vLWY=(-UFArUFp$n+yLaV{7PoYCbPNo}%3T(pxUP3`)a2xy zv9g*e?h&X|ix8Rk;N_Pi^Oj?YzUIXX5+AHvmL3C32tVT9MgkWp(UW$$Oo1>2I3yuy;0GEOWT-`@7N8kZsfg0y_(MXFs|H0;dX?Qq$o0V3=-WN$MFRlUF&o)-anIL?giytk9H#D;{jvw+~fV~5zmBzF(v+)KX&jOaE+BIBozultrUpKz`9DxbnnY%WzmCNe*W+YY|<+-2Q^C-6~g}JRVBvvWSS~No{ZydceTqOT=&eY zh1BH-W;v%q*J(GHdqf7as!R@BrjN~GSuH*gtyl@4;37}ArDv$tNz;p}Dw0J}_>)y- zd~dVsZq_#Dz&P&RsvC8%a7JurR$Jv_$<7l7@#>kvRe=U7$T@V-FP@KKY-#94gg$i` zr@y2*^^#s^lmq{Ba`H{F)b5=F?LL;i5%GGpacAC7<*U?Qy#Z(DXiuay^v6Dag7u_H zxi3CQBwi!KSr6ax?y?JKK^)ya=cr6s8o&Nh5MTWS!~S+?Q$$9(d=VTM`eZl{)WMRg;rNUNZ1%g(Gi-0>Q`%LhJmQK8k@ zvty9Dw7;QSt-Uw>YV|;P9-m;B`RIH(Pw+$Jq{iM#ihr+7wqn@@c%L)|Pxow5-(&0N z3-BAKq1nAQRjP@H^3(vN74C;|#P=YhGI+V&V?QH<0rf6^03}X8KPp4s?x#w(l;lwN zgTZ@qe5^VzSrpZ;iz4HXkGoCvDw3X8uVIWd?=PFIb2*ayMMSi68K%q9O#V?Hp|TJ< z;z?2qpI8?~gFCyC`vgCE*x0@lBbmE>`o%U?<+87zS3IId>AvH5x{UNOp_!&G7jEVX zlf4+BX8maNBRC>PcUTD4HDN{<>bv-jAU%L8z1#Ouu*Dv z>9hZ$H{p%3A|MhRHF#)l+;xdUPjb~W@z_&rv_q-lqPFn*Dw#Yri0XKi@kcO0H`;4M zwT;F_T^o~7mAa|sQknDCQn&LOWN6r(5qW)Vo#V~&0N?esMS7VJqdgl5LvefaFO<<4 z3X~atV8xuJ<3+rkCs8sKKKL|fc&}NNDodpdD9Br~!#KZDHrXccj*z!LR)LEB)%?PwmTMiDPh?}DJQyFmes}j{odm}(IJLLd(I@FK#AlV}f=UaVD|lf! zfzq$y=LDgv{zRjcr^i6ubD>@5*4cBfGJ(Ip2$#N5e0HFEJv-ZLag{+b0jY{IvQxMF z>D%sZe@%xx|B8zYo(q(|;nMW_{oXH%x|Yqk0;7pa8H{pEski&RO^f}#A9|Sph|_?t z4IAyxWvM(ThAHBOnf4uPli%>|ZAVFM{;DApAWuR-{sQ(@WWq${1 z-Qaa%X&HSZ-r}xnX2~zl_X~MXG)%-K@6}f{c~$(J=PYdKJw`-yk`g0wd*@S0qSt=k z)NYUKeLm#;lmpj#B6Ye-*{@%-=e6?+i*SS)m{aC3?s_obt#^!5W0SYjC`%nhlfghw z|4s`_IQ>3%koR^m$IRK7P);~~cHtuS;TI8C^MTLaZgpJ-!rOBlN0v1O=RfR&;|Fcm z-QVBZ;ANpYO!P6~0`?(Z6xdSr?sC4@5{p{hiNVSo`;oP>?7hWT9s4V@-e8%vvO43Y zBSeoe01&0OQ|*6!IjNZCzF#aQ;`Qj-;!f?VfRn-Q*rmj(U;R?e`en6)pmSji2wc|b z*g)MzyF9tdurr>*LN8Tfu=j294Vkyp)aprZsr5VtZoLX;_q;Q}fGI9vJNY0BIvpsHvuam;h=hn?gDAehDTyw|~W`mY)ZX>_sb0A4qKWakzYIYTiGt!r9d z;~VSM+adn~2Nv3QuB{~16cQ}n*_&)lUu8veI>8BKY7GJ#11LH_P#+$v~A!v^|@Hep;0KlgwQZ?&T=M{6N>OaEm1tN zqTf%Xfa)8`Z=Gxx z+s;e0v$DZ=is);eq{u!p9=tW?IQGtE;nCWB%II)Kn_@lLv_r#>5ODkrsDoH;pPQ5Q z!(MnRi#f^KN}eG7DCcO+#o^gms|D?xt+>B5pdcEaO3Hn*rSG`))1Uc7yql8Ag~Gq% zq%kM4ASXb-4CCBMC}vi4V-{gZ)wM@N)Sdxa0dYHZ(1o;0Y={pAhl=?~|ThVDjq zJ3%Eao1Nh~9afmC9_it@AH~A3+VMbk%)#2E?29z0_7@hz=A1-$1^FjhWZ}c4-oij} zkz~cA`8Mio{k$yZeM>ItyBULP_j-!#c4l5eLqO``w_sWpvdy3JvU2JxDn-^1Q;&MF z=Pz4f5?AW#=qX?UiHg1ZwiA<3mhIXJJ85kqX?Q&JDb!yd3{rcZ(w(*fb@zfYUO12) zn5jf#(b8k?C0Tw%M*QyJR?m?0ytiz&QiFTYxt~6i$GrdD__Qce#2d3lDzO;-U+s1I zJa?tqZc|_Hz9r%5CWp|XrFEi>RxgJ`>$K0Jl-EyU6&mpeSAGBTgG1~GX0D}$?Gb#I zsOQyM!bSG1?T^0-$&{BhZN27Rz6|0V$uV3ls&O$_*8dJM-R5QEiAeyrZT@~NR>1P( zKq9dUgV+c4DJ*=4AV?ass1t9!CCfs-2_RlmuX4A%a^aVdpk$TxcgiNgtp0+0dsUd| zKA5)1Q5KyjA%hBbO*#P(j6wUzlz0-zdgaR=+y8dkvq!65i{P>A{MVz z#Ku6b%PpLdY8AV@^;9I=4ub<6NO@WJ!P0wOtqhe~XJlxGvsWq;k5#Of z99pCou%G$5XHx0m4hwm-`I3^n{QI$$QtO$oP04{L5K1ZKJu&IZA=GA2x(QI0Ng*;^ zlDJ3h-?w}3R~c5OjLB&xXNfQ9oc-j3BTy$B4&O3~_B$Zqj_IFv722en`ZrU^NpbxA z0@W3!+lzN>qpKe_*@fHb7^rEMWTuL&cZ&6Z6L@ZS+4*7{-A$MSK1A~7?Srl2_)Ob# zmhL{^1v>QXnE;*`4vLb8L2W6%+Z%2+I=HkFqf3=NG}=A9;H_RyKFv4Zf5?nZ&9yE! z41keaNN?~nG>;GsT`(lW1-}msWEUfR$m3`gX2#zH2;51_gfeI%3PNUkEgYDs&e-N! z_+W8y_vn777!LnR-4d1qhIXyH5~dbJ@ZNg17gz|y%Kxa0@TsuX#)^ILYbh`LJwi;3 z3(V0wH9ER9JfusVFHWor3vWsumpTy zRc-=x0hhy{H|+jcm3VBe^3k_vSeg3{3+q-zQq`7C4t92)UOm0FR>eX^XSr<5E$z*oRMXBq!EzP~ zs>&S8*-kRC=%M}WD`Qfz9*z_56FHcJ=8%Q$Ds+u#bK2w`DC_FSvodd|fS*cQt$fUv zqiB{Ix!YIr)^Z3b%d`sTeKm9=Yz8~FycvD>jP@+R4LXDKSJX9}gHj`oz;uNM84Ud#z?GKl)V z-o819lG?@=F-_aiO(^afq$GeXeSsCyk9dE3_Tu^VeW|qtamj_UkNX>D3sKEc4HBhl zZpHG^EARZ09;)lGpI7&-sri=n7x(!)CKJTvS&Sv=Sqn4o-d<+SdibnS!}f2B!DO4k z_bd_*Swh`gnH<`7_BN6~^UEs|Lgaa~rxGa`mUq3Rmfy4$ZxDXGX)}onONwYRrrZU% zCn|v+o;1l73)NXEI&wL3N7pr{l699dxG0?5kBpn&+Z_G$akhTzu2e4z{Y=jbU)SVD z54EHgy(l}KNGpZ=Bo~0PtIIJSy{#FolUpI|>H46|k>rvAWg#UuN>CNBr$}dOTIGC zw3PTop02sL&8SX%Rw!wy^!DV-tyQl0UhpfcTz#`} ztUTq@f=K%Wt-u%ljXf@E`cB#yd0kt++H$O!*}2yb@1z1c+<=bFJm0^Vb*?e*1o>&hW@kL z!k3GqOw9;UzEG#w6R?;=dhDi|mQGPlPEi`70ax#{Z2Nutu1MkNIv87Z%Z|3kmnPn~ zgvZYKwMF19zJxYdK2PRJ%poPLub!aT5Rcf<&nSf$$GQnpfqn_%a&dV3>YKRF5{czl zN>g0mS3&beo8%?&#Mh^$rp^i|pt848`e3CywH(s~_;`DpzXp?E?9lrD(f-MZ=i8=> zF9mF9`ABnJc!EoUQ*U#{ys>$<&cM~oN2TLPxrMl<8ZE5cI3c*n)P!JWLJ-mjaKx4E zCOpv}J(BV?^8Cz_+mRGb>;|(o|9J!M_IjBU%iK)yZ z;Y;z=`?-|6vMgJEaboS4-dwSKO}Z2?;mT}L?(gd=mw!i~Llw}6M+6L_WaL(OPy5zZ z-n*5JD8ZCP>4+)A2*j=WCnQFNXaeJUeW%IrE=9s{eK%9cupp(p@U2;?QbvUVqVC3V zSXGz8t+FgueV366Isy+i@chDg+I#G)7uLv+$yFxRUK##IQre&oG z!S6NVYS-e}WV7Si&0tNDtjH`;a%|b`Jep@Bzni&L_#C;o>lX)V{yZ{j$P~RnWx$$O zJW^OZaz^StpGA|+MbfD7A*KBM752Z}^`Vv|z!xzY3xCiEXr9Z@iVM3$|MWt~eNWj= zY}XOzwohEDenV0-oX$aG3%{zj?q)lFxbqb}q{$}Px3rEVCbh=2mA4)7(PSpj?(5{+ zU#Wkbh09LwHTRnbC;ZONEg%-xkAum^&oA2ZyxIsiWGbufs;us!mkPD&W^%^0wg?_f zaTK=DV(X%lb=bV};l1G3Cn)+Q$eh~}Y}?u+ zxxN6dc&4y;F}JGr>LHmVj`p?klzjgMWjuD^9^JntPk|2V-G~?07q!%I@OigFZOlMFy z`%;cQHwW-DWnb)~PJNG&ZlrH~wC~s6<*PBbDnGnU=N$t`)tHPiGk)J-a`h_2egj1D zzZcF8`98|?)*|t{O?@xDO9(7t*X~uYK9LCKJiv7@MyF|3P88K7E!o+&a`I2#E5I8S`17 z2>Q{H$-njMeI5@1nNzs3-?d>7AKaM(;LurUlO45{GL?UWv4XQcQEZAKPyfhB=ACpn^dG7p zags7vU>bMyF_T)ONj`4@182^Nx+UuAP|eg#xykCj+lt(_F;ofBA$^CXC7KOesHAmf6iRrzNgZ*NA6NN0x{S% z@rxy#;@Z2;?k{6=CodS0+|%%Jy;TJl+tAynjTFEbN&pC=uGrH*AuAuG%qaIR_rB+oZ83>$Q9a~Ou@SB@i*iOtNiJuJM=j`r$;r-|__nC)^K^O( z3)O#MVmM61qYc4FWR@7*G&w;tH9;e!IiYBFmiQ^cvYa5Ft&tt0Y@#?;BhV5UJsEOI z$%(!8hU8@)O=>N8UY(?<1UO9|n2|`5&fxq0D)`O2x(hzA!RJrj0{y(aRNv8TA)H5uJQY9{7}&wet!Meca*) zJKhFse7yA9m*935NQJ7#hf8Xx>nGCVXwBnraOX@%W(00WJtD@F#8Y zC*{n~TX`r5s`&RE@(lrIm8W+-k=RtR@rjJ2ST8Eqd6V*zo1R3jNJy+GUMIfpy$1LY zHMI^EHA1@@c{`cQ*1#XcKmzH`fJVPv*}vInd!0(wih|<2r+@iTt{3Mdr-sLp_`%CI zP3N!nohN#n3Y9ofm{+jSq*cTbX{U+_uj zT9cYv+At6v-mMw8CFESbXk&3}>)YYwR~HmDss`SUZ<7-}e(>|J<3y2Xms$or92Qup z0_N;G_<1UdN5FbzjGhsdQLyaU5@pW%|a%!QS~U+6fZKPq9v#dhI5;TMF;BO}1c z^q$gmCu0YH$2XSpAtn}Z8Eti5Pz<~k8fylrnHG_XAeR)ZiMVUC?$uoEJyx0lNN6LI z!(9vDk(QF4%d#olD_gl$eX5))YOz_mu(XV?AWmi4#|%x5S3dz@lukq5H;nE|Kl3Ne z?%%X`iLp{GSA6lz^GOP)O4rOvp{Qq@V81U`q15{CxCaNcE9qKzsmL*rL|9c9DHtkP z2h9D)%YmO4X0}SqXAKJW0Z`fVoZ_YpteX*&0N#7uxrhX(q%-kKrU}@aNo}~vHNpw5 zzlfBO11i!A9EHhYI$e^*tOl5qz)ziAt%hO?ZCWL4vJxwwJP)Oi*!!u=vNa~!!A(z< zG8;~v!p)&#aK|{6J;yvOLyv?g(hWHv5^%H8By1Dvt*!_zZYr4*$7)ucRrUSSvqmZ_ z2d-EyCGeN55Vk@-CFSRx==p0>!=yLtQm?CsTwVBKqGW^$5q;+~f1R$7nXYi1nW3oM zqkKpYVg&>#3+#BJ7_{;TuK0{utc>w=Pt%LbP5xaGAPQ$*WMSG5anLIR9xVM)gsos2 zz(X}+>bcqz#Uy*Jp-J&xvC}UrrAizPSxe^47mE@rzjcSxG`;BL=9V|O8!H(UW&t@A z`C8Ir3Qoxq$EFweaQGpHG#jU$B9Wy1c?LJ{iB9v&&Lus z%qTaf&NwK4kX2B(W@za1iOoh2Xm=)ze9yLUwzhNb1_UiCKRj+;@xGqHYx&VwO#SsD z35Y_?FdfbCWI)hxD(!BPJIRNLlvvSOqFUaMs2dHb8&S%fAin5Czc45`1$YBB>Bx80 z%L5a?q=pP5iX`~WUn`Pwc7E-gsC7IKXmc1@g(9nRrk1{`zW&-?BmZsQ;hxFizHyZ~ zak>-ueYd>*Y?kXqg^^)N)fHdLDp)((C{~s#8^PeGI{QFP2TO;8$Y7Q z{)i!@AF@v5ZI;%-#f|ubW67?VE}fWe-J*Dgg>yH**q2ljMQ&~tgIo71E0{c;!+kWs zucuxjHhh26(93($h>-ye|=d6FNVv&faO>3*v+~R3OjV9N^!;7nIPa|*?@QO@@qhmGRqW-FMSk)d ze$Ou|LZud}<^(msRZu@;;rrAx5E$9nR2^D#H`@MPW3CURs$Vm`o(~Mb_}kY`kNi7> z?e{cN3&q>s8yd3_DUl-Dx*eg%tQ%~l7i=}0_f8A(9;;}2ftFSvp1UbQm**X~Bo5*}mg9NZb)gK9=E1&w-9ia&`u6s1=BK zE((e{7YL8qJfvG#lcNl;7VD65eCIv=Tw-44+2tp05V=m;b&Z`XX&(6spG39N+%&8~iGVUtRsga=tuwB)VI4n!VQr17(^ty6B>6j;*vvBAs5Qd;}E-sY`F|YX1 zu$K^;l1y)Uj$hIL@@-)~9phVk^H?ITD+uzI`~dFwa@j!;3H(Pp`dS;RoK=d1XA3@mlHSOq z+_if%4d=Ti^3o1Fc(im{AL$}L^zS!A@^RmLZXW60k(cANpwd0)!Pn#&5lzEe{kE|7J@_*gtX4{rUq0q)pD#^f8>24XMJy@b}nn;Cv0~H40enRSKN_=)2th z1lz@OQ~;NlW1=H_51?kgXIu0DDS@J(3M%Z5cJ8o}p}_HwwIhd>dBVP$ ztuE+I6{}OjLBhSC3xZPs2+IBr?p=?ij0TNThV~sRcL{Fgh;i@qD}g(L0T`l;1gtcW z7-c#FpiHDB%bdBCrmd~n(s)W{LndUzCd(y;4wMO%X^<2hw;&Csj~{_(8BCbWU46sw zWcn;tqheViHXkHHlSSSr-TM~f=rofMD}b*nI=tg>+5gArFjQVV9cHORz{{s1Z=PF z3N0BL&mH}Azw$}YqAj%yaHLC+Utl8OSB<;{3gM3r^C%>`iPUHiT}Wea85U`? zlM{ME+}ANDG1r0TSnAM0po;0lJ};-=k_iFbEDDO?vAj`{5%=|Oy}BEg7-JUlzq&i~ zXsY-B-|rGhVV4Xc>;_ZDj2WVB3{hrc6EbF&c{plI=E#tQgODU4nTIkZAsd-z(#}k5 zv)}9Ed++_{{&Cm3>#lY0x@YaRoplcT^Vy%z`|}!}ulMWm(ld_J_d!X@%`13k^oOGB z$t~Ei=Bp0OOo_fUDvq|oYM*Y@K-^qe|NX1BIj>rK0ObihW4Ju!_|emp!Fky&EZ||l z7lLJ*lN!!Y^%%B*1fcSVn~&>u_?fwx%ft%hpM2z;(nm2J$3dDUhsZ!DN$~6dMCzmS zn?oz3Bj=rxh@?6NRfG5Zk4HH0@=8H}d8fxO_=yi9D%`8u zcao{H{3sK9WX=DhuAco&@QRw5etck`3LgDsaiy=cDqEHx&gN*Yv-jp8q$|$0N zQqvk``6cnVv3tZ);_=Fsi=yQ?iP}??7K!`6zY*oR-7Y0>JwPnbH#CTk^k5s5-G^s*5Gnc0j?$3KIameHg*DT7`oQaTf}IC%Q>OnKWKVGy z;KH=sKj+IJK<>hMcyW2^tqFk^Hz3uWrovq4>W^hB?C~kDYJd2;g(a-E?p87^xvt6n zuE}RK{8GKhXCPV|yxcRCl3pE02^~Jns2rvclTSjOAB8&qYPE&q(F&hcO+y{LQvS9K z{E=~V$!lZuo$?-IX;IM^Xcx9;U^Il>9dxH&@m2`npC@uI6nQBW8j}bRG`7|%$f1f{9{U0Bp7Rr$NmSAlh?+kS@ zgI_7XtXXuae9ojDslTiM(HlyY9?V5Tb~>afc@Th{j%tb~QTHASQPT=B2Os@^uZDn( z{DTD$75PKk z(9BAVuR3VUQwx4;Jpxa$PcuD6h95LMlIe!`qZckbdmt(Yw&oRW0=D*W{hgKTmOCr+ zV+QIvNNLBM#JX^vRVIL52cwwM-S_KnjpzK{Q$zEg9qq#K`rale>`!S)`}jO`X2Zgg%g&hU zZzYYsJnQ@Ce*I#ywqbz++N>1|R9*82R*se*J2Zvmm#nVdNsQG`jNJf(NJw@O^(+qh z?NUZkjUgfLo#QLDpXI$h6g{@Fpco8_+w4YS1+HlJwSVvF{C@R&2LXBi-ouzXch-$T zj9h>cF>5F~MIQEn-;>=Jw4~_5nqRUus8rT0t<$hAoR7zOx117B^m5{Y6PvD*|PEOv2*BqTDNMAl3G{Wjq_}ZzrPi1Q$-XQ!-rshvxAt;{~%o_D08<~ z35T`z!A%oIj3NiC?7(5n?madhR-18Hv3oLa$sk6Lg#w;uQsKzV&fO{bNuL63+{;+~ zxHNslM*a!stIYH;4HM(;z;vpey-EA4ZZ?BlGCeEQYH zCCvZvH>5D+)~3NqS870-6qvYp{u`kfP03G`XMW{&UA~p0c19u<;AVEA8^q`g3 z`>5+&`366DVPbNUPJGf}WDV`NmsQ%>#bU;BK`gbV|l*LekY zZ%@?>zRXeU%Oi9v;|#ldh*Akkxw<@?*%*N^J2uMnz->E>qNnGEMD0pc9j{J8HgU zUcb0PqM@jiJ(BCZ%L9z3Q{=~=JmT7JVNZAoBH0I9J%-e8nL5d(YR~oUe!k_CN+9+) zJXXF(;O3V~xRC|w8bv9zz1?B;{{ZJT<~>D?X{>xWr~9^R;U4tv$H;EH7+XT1TF~~C z8aW7N-CQlMRT1uiWly9rfm47f42mXDfq)A$cKd^~P`}(+aLj*c$)|!So;;e}4}*2) zB6M|8lf@>(;0mFKDofdpz}v0+O%#r+XRjA7n(*bFj(xn_V`IA>Zau1GrCe>EX^PSd zw+eo6({mPRW1bs%@Fr`#rsK;kC9+c1;GAu2a?uvh`DIxM4QTqlE2Dl(hYrW(bRJa2 zB$O~6k_Dp75!9zbSQ|bEPmJd=VsZt4GtpcNRnRp1D$E*>EBDV>_V%;${})Q8+)PB4 zU39CviCgmbcY0J9Lbon-80cS4I`LRJPy1^xRWh#NEY{FNE%_Nym#RQb0RQYlwh5=m zTs)RC+r+v@Rd8B43R~&u`_kAqHb=YBz){R0!L34r5<~Y$v*UzND%m|t%dEed<>lV< zde2*}iVRHz4630A2qkoH`l?aljnwD=&>g0u)*K2GXKvMm`Ok{;d2(rs1wIUd`p*%rQ1&773ci^^dq()hSJ_u17Q zx#6eR$A|OqU*jV9YS6w!)ys{?`S8L_k(=1AiQ7F3v^!t!m~sC(f!~`j&s`= z?~KvDM<;Gya_8QVK%v;2kQo8uOpz4ui;_}jBk@_nWpZbtnq-1iMLA`Hm+>-V$C+em zm~93PSC6)Te6#QRULdRGVpj61@uJ9ae@9gc% zb3A*eb=;-s*wJ7^zF(vroj-HOUsxZ_jWqp)QoWFkRt^#3RgBih;6-oixU-Q?s8}jw z*orc*QDk8$wK&euS#Y1FR>}Nq`<*eI3NsuG`z0B}DJ}|Q@DfZ>%@(n?4}_TLv3ScO z3{go`{4D$q*fB1l;-VH8Ze}Sl&YyHoRa4J*y;#Q8cd*Z1VWvE${U>)0R&T<$*xU%k z=7S(@?z3ld#YvJU?e4!6iQT3xDJkj94QFhgU2CjOgYQAAQ3ruV4c~v@-Ovs<%-Is{6m*FpKdn)<>`HvyGF|n z;6p-0T3T9>fo0JP(Uq7*DEFdJd3I9AzYqx8rN;btTYGzEztWNt5vEgA7yy%qbcun3 zME!Hq9F}4T_RgI<7<_-JLtlSB;Qj~^(93zEfC2ds4Vvz2XetaqJUH;)a|LA-6coy) zzQ2-f@Mt9^C2d;>PEh~oinNnWER>DC*Ux|adAx=+n&Mpx_yTCYqr(pS%%Nf} z6N|v1%F4=+p3gdju|TRCYwGNLAv~3xgi`}by2S|cyBuLTK^*aIGYFwU!g28L=2Q$2 zYx-x`03>M|w6_jOpO)i`NG5|qb@%_J#mB)+1*hS1;Q`n%id+t>mz@~oE1_N;LgXN} zKYJL*{lE7!{^|hz3r?%6M#cHa`J&iqKzD~HpkYoIB|5Rg$Xo@q1BpX87GykVP^uM2!8pK}XG9e8U0E!w~@O<9-5!gTP z-o3NxqpVj_AUn<(hOl!aB)r!D^x6piTJZ)70W4>@B#ivV(9qj~+yC8_1?{|4Ar?R1 z86iZY!dg_b>-5)Hvg0)G?v@G2%I;l&J;Nny@Oupn4e$WKy1JS?e__-g+ayH3H*X`? z*3l7tprXs0fXpmdSh(d`gNoe9V&CJQ`i{Weikvfo$(}cMq)LaNfs7-mo}7`90R&MI z;Uh0Ka-d(8i|JJK-F+Zryp(XpO=?-FCkYhf=U0hpp*Y~$0Y#Jx-w)iba)XPpsyN)J zQSG}j$06}Q{LrHe9mnb#|G-HBf#|@C$UYOt0mBE(4z~NI4rbfFe~)idY2HhWk6*rc zQD^{xi`r}@m?$YVy;&7?F<*T)xMXxRV{FOn(*vXGc23V<&hS{A35(=C05V+@OKF#F zk9{e@AxgO3f z6m8IE!wH|9h)7J9&i7A_L(qV|uTEJ6d@5vQ#Yski`K6^GV0s;_n;n#9XV0jcZ4=U= zHPNc?bo~y@YrxkHg_6G1ROtE*%Z72N-J!Vn}QRHv|t>1;K1#pF%Ex zu;I40wz#y?b~noE_GaY{wg5i63w+)E(!v0KTsgQ{UcUa+N^_QsR>lRQd`_C7>pHaef`WpKjbEZI z_O{o6(2mKVA_5z5Yhot~b4)s|}JGeb{wLWoLd5RT#S1M!m>CqUjrEW;ibJ>}fEpFj19 z60u+u&X;~|^r~Guk5S}f@Y@)$d9HH}W$7)GHad#KG`=$`x&A?jnpL%Oad8o@J_+vl z;~uxpMX=1wY-|iznZG|2R)YVY>r|vcyE&Ioqj7t zy#0kRv;t^7K?(eD4qQt?)&W7#?D~ffJqq3eV?C3)3LPCC`^^%yvEG5;#UzGZ9-}AW zIICO!vpU3S3|jLt+=eHAC!tZE9SEVtkV;m@Ngory?x0;p6j>icI z3NEc2vvID_PT}X{;}aBQe`p66mr!5gog(YAYNf6o!BiN05+m+3%#8<-T1))SwKFeX z0qQW?TOm^(#ziVc80D4CU7>@YkzO-V2>{(TX620`!C7V4g{e915Tt;-&a>xzkSQ=`nMeY zWp3Tm(=%K|!iw0002UFXQvC5`rek07UgzcAq`U*`cCmMI61FZhOn_X?5OE*tKw&OJ z5Ii5HJOvP{Eg-~*tq|fat0Re2WBRV`czw@>QlmV_Xgw2?>Y#&tAT2b0vEgNx8|d#> z;6=EQjeY7J!B&Pz(cp||nslp}Y|5c} z5~+EaEk=dXH7A{XYIekl)YyBXKKW*y$NZl^+knq~85_G&1dKPh|IaUsrOEmMYMD;w z0@w`WYFy#wzv3*;AGFEjbbOd4pYb_5I(|04O}Vn#Cg@i)8gh~mOiH?GrAMBAKJ~Gq zz&Kn2$`nibQ}2x@e9fY}d5`2J6I7q~b$;Z0=x;C#hC3{u?jMooM>DkZ{Qq{T?Ck8o zBb5Nq!J^XAb-+dT&aG$M{vBbqYvKP54X1$J4Ga=xrf2W5o}$*bALmWrb=^w76G9IE z%VkBDOaHye0D<5(ULTCiqJj7c?hv58veA=(W7_R447{8W%Umou`$SF;O__TpuR#jT zv!B>ms+)pm82ibeDI%an1SJaJCkPPLsxk*5KrsSFC?i%>T>Ls9lL$ANC51v7L!YJA zv^TTx!GT`^+Y}K^0F8Ix4YwSs?31%6Po6wjO{$xoojr+JtRdlu)1i)=I^^zf)Eb;1 zetBna3Kx|w!E}nwfD2{w2{?5krrUoE>b7f^0y<~h7RKt#q`=|kptH5%T0%%TF!Z4M zvl_JwbJP=TMo|PEp+4Bk2to6apKPwN5@JHY7R~{4cVRDVYHBLjZ*2tz_P_UKXPB?b z;&sA4iWzBXfzPX|stN&xl<(q2OKd**c!Lga0)RVDV*#sq`|yB??sg?$*TX_Ob#c6} z4B6_f&}P#&;X(CNh;FhVFy!dG$I!vIuzMW8v9Yn)AEyX@6Z`bbmoFbbehhnm?jKsM za8HQxrCtY<=|sSNY732q6-*>uzi=<2#I;R5gszS3#i_|Ud_wIzh&p@iOC=Q)7U$>V zf958;Q76L!(KRqQiGj@q{*gwRgM-ExMQ9sUG&ouBJqiPx5C;YEu+5;|Hr!n7(@j-+ zxFMygj-J3iv9zIECy~_WI$DMv`1n+a?%0CauynpGwj3E5DdK2xWc4sjpLT#2($%@_z`P(#Q-GGjGKG(j=RnIXxm5r! z>jtI~$t7U8)#et=lv1F)V_-l6KqA-ZESKC?*oadMBT6XN2|PAN5eAfjZ%e(C6nL;R zrhaY-VGvGoflczXyOFm^$ch6D3ZO(ZDP6&k;#SqnJcKoPV3iUIxv5T`3<6pqkZvK6 z=mxH8(!Rq!JO^r?MkMso z20CdJ+5$R~ZFraeth};q_VB!E=j^*g{(yTNpnwlLjK%s6m zp)fu%G5W3~)2aP!aA;HP2mOHixy=pXMhSp?VIH3Aw*@O|4SE>hrzNw3M^C~X%51lm z)+qv%&mIJwj`sF{78hYh{}n;OF@T+c9}UP^6|CKy1S)gj=JLJtP6J(C54ba9t}W+v zoQ`AvN3MBT+(pTUACkZSFZ|TztO5 z2`j{KBW;3puhP@Uy1HWLEBz%Q90im<36HL*s6be<#rb4VUi%AGhMbCqP}(}+8{3y| zJDBYl7tdY5quF9+!8_7p!CmoRyLOGx3JU|lrZQr`r^}^7iunTBU!LGp^?Iouq$YA5;I>apg7Qh~dhFvSpe_bwE>L$FAY9D+t**$}}_-hUDrn(jRrCNb)rRZ`Mzb+sJ` z!hplg*Dqjxm=BaFNGg+(l7gl$vfDX21^|(dv@~`w10yw%Y$9m@;Gamu3JCD`?**Q> z>KYoB1NsS8%bfy!b#ln3ff?)?+JHHKeeJsl3=eoI7laJ#$Zi-N9R)rkTu=X=Tsgo@ z0!|C7qoc!pws_zHl9MB+Ik&#ucojkkAd8>IB6Ih$aD_o2#stFHhBQf6Y=M4fn(WTT zkTWE21O0WGkfH4Xr%)QgST8THz`*_afs(QZ1KQf!ZfrqAP>tcCK-GH`-D{x0Cb!!Q#%xw5*tAW*U`$B=jwCJr{n z;R+2Zdw%}(1bZ4C&<3$ylHqk$8aSrU1L0C*Bvc(^oRFh{{f)>LJt0Ruef`m(1GoW8 z2Vg~C>*p_3EnO$R_9z9^R8u(p^>nS=ZUp=bga)tz04H|4b*ml7t)Jiym{xfLoT`z9 zhMGFyk4`$^o7J_n(q;WTWi*h_tUK6S3PHozh<(WOQj(HpJ*39KPb1nCnD}&<9|DQq zMZ8c-JCK~~+Ge&_VuJYakON4=utx~TJ9H1sh06W6Rw3!7!XR);WH=`T{h(+epN-J8 zVHPKN_o~cHi=HYZ)IwY@_=+)@GU@8-YHB)y(9?p0VMdGhI1K5G@{$_`dM%Q~bKwnP zsjT;zrN_p0_4jZ7&Waga*l_lrVgU#7>?kj6@>ch}Kn~@mbG-`Ou?!AfbL&}DW)FMy zuumr^Cn16r47|!5bsZ8gg26PBiEsMkh2GBAwksV()>f?Qr4@bf%-E;YEaHk;U5#r5 z%`;fa-0pL71H;37slz2US{4&UvUtR(^Cav789mOh7=~V6yr{6o(M?TFNXpco zPtX=FC@9FX-fNv%ElN3OZEd{)gG@=Jn3{xcVmz+VqKNix1FVG&66qwSc_eoq!Z1~1)`R1|iFn_KmrvbA;2-${VFQ-wB!T}Bqm0W7X~Pqkwx4L{h3hkBh! z$E8V$adUS^g7L_jL156;z|a}Pd-v|``fF-N4?MWOlL~&^&u<5?{G^2;(sBcBx)=3S zYB=aBtbV=bLhY+5xjmJ zrGd0=%JariXmzB^15kp5l3KZM4IKntQdqIedqo>!Dlt}#jTXJa#eg`y) zFhwaTspb)$+U#tlS3{_-5`6Xl^pqgKrd7}_Y6Xn-C;|44~An(3vD|kJ!%fM&249gx^fLsEigMEEoqhj$#Gtd_E zQkjI*Y@OWqVZh9SrR$vPfWDm&aHy(E zhlylnULFxdc^1F*HW1&yd?7=C(%K3@4-P@YG>fE-ZS}vv34Y~84OQdU(VyA(mqRTJ z_h*ClZG@OWSD%&mu!eB{VVXe_pV8G%wC+SgHX$P7DjlpGq_GmkgF#c`u^F2aVtpnT z?|(r?SV?=%tMYBZEGQU+s)K!sr28g|e2$joFa;Ts=*Bffh6A;WqlT?4-^ImcnXgQb zH-X^!2*}BRDz1$yzt)Ig`T!lSNV!4_Baz@#<&xB*wZX~MY?cB2sBx{xP@@u*xu;yv zV<6fL6ZYywZAxL2<^=FER%cE!@Rjdv3}uJl0wE77#MKNsmBITq5P*%#UG9l8kJJWK zi7Kwq!IWag+T{@Uii6qfqzahsT3lNC5WXX%*6J|T$iZO_@hkAyK-!#`$Pm-A26_4V z6BZ{>(E`zhG3ZGM2TWlpgwtOEUwm~j)wz<0Ov5GheG(a7v|bowdiSm$tYlDyAwmGy^)2I1f5Fs>(@`GbxZDVx$^wTq z84!`QWQgm&voOKq;6ENbiJ65ZiT|$npOib<+K%wtcza8(aDt3SkL81idAM^A!2h|r zX&ch|$r_))xz~-ckA>d`^WTO{Ls*GJ!tCtl+Jz%9$~M(|R-WSnJcW%uVR6VrB$#Z{ zW2u&eAl^mx+7TIfMFhe{z}2(uJicLUJP3aTf>7tS1T_lk6GdyL@0np^@2~OCkIBf$ zn4^AyI-w=AMA`w}`sctvazesG2?_At_f!hhNz)g9DSRY^mN4-q%GNx>^uSRT%T za#)kKmsxxnDGXch!*;@uP-Sv5eI?j$||h24bb5t zl`X+)fZBjZ!!8M2VcDBEw)460#8t!06Hxp!Zhy5dF-0r~Sot^cOks+n0VGq5B3PSh zI=7xN4JRij3vtkY_m+ajiY!(9{Hj`}ctAD;AoMhsTykuz@cx(}#F&RWYanp9xC4Yr z*kgwDT>gLb^Pp}|3fUHcLQgx(e!YKM2a9Y9yo?l{{@4F{1XZqNlT(xn`+*>MSJG5` Ik2QbvU-B=fL;wH) literal 0 HcmV?d00001 diff --git a/dev/3-parameter-estimation/ffa91456.png b/dev/3-parameter-estimation/ffa91456.png deleted file mode 100644 index 4dc17e533d368b2cb5754804dcf08501f159df78..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 34990 zcmd?RWmr{T^e(zk2~l99lr*RyARyf!tso%{QYziuB?cfR0wN94jdX_~NOyO4ci*`- zzw>|2{c@lC`5gR!%Gzs=Ip&CWykkuK{qZoA>W{1=E48-sbVlnN0eMiD`X9TFR7H?p`oF63an{*Blf7wO#06vGfPXf z4g^p{+JFfIJ^j3m^YiCk!|IS`oNO8cGjl|X-iHq#-WlTFDjXO8hc5!ZG+z<#|ElkZ zTf#%i-$ol8S!)hIrlK1q=3V^m&Gd)KgGzjhxH4Di!B@)fKv z)t?Wt`r2yS(z8tmiF|4aI5>+;ZvK@qkhs>~;p0Mu0~wQ(sHBC*SUH5%cnd3Plj$pr zQcqdrF;KuuhCk1K8B{jb&+2UrR^`y%S4j8i3q$!3dUoczti&4mD$MCu0s=nse0}?l zsNx6s;K3vO0g}6fRr3qK%F;aRz|##4{EzPmr}c$LWQl#H`_Csx`bTrmSh7x|iwZi@ zWJ}lf$#iU9|Fc-zaMTf*4rH0cW~~%-{&2k$hvoe*8|B4Q$>HYQDg)_*4sI>NG}Tm1 zYc-z|#$329C(iu$q--(*!OcPS*leTc>zBd9AZDI_$y`lNax!^|{I$O> z&+s#{7@3I32bB6h88N$ky^R#3zSS$oU#Dwf1u|(SOuRCaeIqk*mnePryuW;QmU0rw z54{%7w_O);NHFLbt8uC&u@HHO0#~Uuh3oGb|FdUKA~C@VO!y3mIlqvHj*zi1-)rqJ z$WvtU5R<3I1uv47#;9<7?&TBuQEkVJ{Cer~fzl#R*P`>a{(_t0ghJFr3k$43LXCI2 zdN|~`>o(W>FzC(gv*cASq)u7)q|IdJ74L3nnVCdEIif|wy0eGJ5Jz^s=(0_QzM*Pz zSgtIy&_wbIx9Y#-`?K!kL`3W?lM#})I)afllaQs9KB|)XBj*KOTc%0Xataj{dEdR8 zX?HQ!Xa3N6|4Kw|U|BC`qXMfrQzbdfqAh1&1&d4vbH^&`2V%djE-a-LgriF5a^RM0xr0+F!)we?x$V6MFZZv)FRrdLkyV`PZ2 z@@6*g3)edy?l`{w*Z9>(NCau8R^h}0qg!nB5ouwYOzCWzR)*<_rMXkEuCbqR{Y0#o zluug{QpjV|eYA$W6MG#pCN0nG0`;V| zuXz84yJk5@zr(fVVVq&E)$%G-96F-6l2~HJvOMAb@J7@e7rD259sJA^5b1Qt=sI$Z zSGglzKAwB;PEJbdCe3VtRjiNPUEQXU-D!M5W-Y&F1y7|!TJ7{;x0BA$gOA=@a(8zf zZSfe$+rJ=6U6}ON@|=*Z)%HP;wt}Bnwo^XhdLIty#itZO zsuxREEzCT{gk+C?63t&uPk%eRcS9hCwDBvtv?x-wX7#_N1_*tmYF;wrRa%w)1jm$d z`vV)g&;A{&`wS#c39c7?m&v)Yy1%RMBNe1s_9j4yx}0Uub#Z$7jq3^B0_h`P>1#_H zY(~01GoR+JmV~&!#L!NPmE00KI(ioo9_jeROCTS?wR}-ZT$#WJ8TmB7YxVV8TZ}U< zCQ?x)CAs9kxtoi)Ui1cU^;_w~6|krW@9T-&%4x{Y1H}Rd>vm)mVt%FFg5}us$tMu8 zy~G(Tj4pPqS_(uT>l=>oc21Cyr9QgF&(5u=7;P*FJ^-LX;uo`pX49I2D1(hf3x1r3keg2M*bGx$2)HN%0Qx){B95fhxH zHlq9fzg6W;4}=Vq83oh|1di;Q8N3zz=1-0CB)^MTY;7!2&Ox{BXplxuSpZCh5M3HX z6ON)1#~2?c?haB?R@-OjFA}n@?Zlcrqj-uu=lEQ=1;pSWTlC2NA8|D}w4Vd~4;I0T zuJNC)g|WV21o>0#SnXVsK8F}d*K&ULM|!6Ty0sQkV4u0yFJm)SGF43p_w>l;E?zfl zE>~cx&dR%L5Ec3Uu)2&I!}dB1YB3eBulL^ko7LME-A|}_FJQ+Zl80MwaT)<9xe?ti zS@>Nsqw}1ITm|Pb7o&L>N*AJx`D*(A7m<6nM6uo7!3tUgp`M=|_J`Ire&Zk70MZ$d| zE^MuEl}V?W$Mu4j7GRl8xx1^KkA>C9k*c&h7%hVR^9T9Uksr1dbDDG!lcQ~0A<5r< zvrK}OnEa8>sTnw@92mdBtP!c}sLZ2_P6h!MQ?Ledl zNRhS;N9|9q$n!L85}YKpbLy=;a=*OZAx5n&I8BT^7%yoCfxZf$r2^5E z@XV5+Jzw%putP&^>_cm1pvz@Eq!OlYxO{wTdblEYM0n*^pphozDI-r4K0e-e)RUOV z)#M^wPJr7cJj!{zbKCBwPBh8H9s#+nJk={olmvY_sR8*ZcT$`Ri2f z5Z^ssJ!5cHM$}J?h-BkpI6wc5%SmgB=Dkl?NQr+luz{8XP$!uA7Yk4eO9T%`j1Nj5*rh zI2_LJK0`+<=d;P?U&ls90!|Cdw--LHm{+b4MdSjtW=$boBy_1dL~J!P%6b? zzx6U#0H&H8=21|)tY@#-?3_?J_wZ#R(D3;Lx9E+-Af#!$v8~G7nL5K%P5r@@(Pc^W zFa%>r(pf-|e|=EfBZdlv2pKDMA6_5o=~0b8SAeynos+I^_gq86teYkfUhUeiv$d6z zTvk5JK1Tn_pkz+AbkTnCnHp*C5~jgtAkno)I$Az<)7&jC(yabftkgu{{AoIMRrCKq z=%hBdnT?0kC_&EaoD1HR^*Q&?gk-ctBMl3`-zKv8d~)LE@eRkycR(7VT6#=h$JM0j z@j!@z^4Q7Tt(1E$fxjo^0x>8kmKejWVkoU-t`~d?EEYs+oFhDo4#6C*m|DL6iJrGj z3DdL^pTu3ueTi`LA78%`Q#txD>tPDT7Gryk`6f58tc@1jw*aD9$53{c9|wU_!0Wic zG)iz%wdnkXanh-GrKsh~!rAA<5JXU)|0CI-r zST+k5Q0OuDR{6t+7QSz{Bn!Kmsqy=||4IRCbj~AR_-PW~8w}~~zz^=}BCV%y^ z<&Z74ga|AK1eE2J=7o)kqO>8MwT`1D*Z5RdZsu!PuR&lZJHYDz&~-k*?GfshX*7c_ zuhF5Z@wFQcQDpV6bWuEisW{-Tc=$;k28$sg|Q+vD!k>leqUIM%^?0K7{kc0;dNm$v;qohH$znqjggU{y$|&h z@elB=9cv+Syjx}@uz%K(gHU|Mqq)18>;3(UN}oJ^h$h7bhoRUpW3?O5^nV6f!67zZojBIR zUgMmDgN1R+sq+_rvZ*Cu0r(*vc8Mv~@M zIZ&S00x}~FhXGpzmV17#YuzmH9@M!YZAH3jiKY95`i?+q#ITQh1^&xcx-d1oj7q}a z;H@?$nqK5Ss@Mj3oJW)VKl9F>45}La%ePmmH)J=?mRhP)>f90nY?oe7Vdoc0z7g`+ zO;h?v2K7~#E;)o!_FpxhEI-Lxl8_7v*6$*MB9LWQzMx6YTMecy2MAI_b=;&;%MU&r`Yr`)5m{rgkm!1)Cc{dZwBWpge#r zGV>|*X7@yX5zxewyNly$<5X6fB?7$ZIt17ClI59b7CqehwIM1h5Oc{%KRZK}$b=aQ z@1Fy?|AkRS`C2QLl8bas!H=qPR?K6vk2Oi^!Xy-u2ohWodl#E z=0r+s1iUzYli9`7ab~2%p|5*PdFMDkzkEm>&M81S0?3p1r7SIR?g^^OQ*7%W-J^Xe zwr}OZ$z#)$pi7G``yJB7MxI!q6tmD-7n$=KD}(*%p5#kVjJRH&rdOU>f4>D4BK_}Y z%jBQa&wRTcg3k4c-VF@`@f2kCwsJbRWe{E?NJFgnQBP^cufBJ9IPU&=2?syy#i!fW zg3LjF&Wm>f({EuNnCL%i!PIMc5QabJ?)oS7~XN3A$*!OZfPT9>ZHD>tc=>f@YD=o^yyF zFm7!ysktPS$z*Y1xt`g(TG;iw-~M1jhu(fOpBSkJ2K#~v)A<3doNO0 zdP)y>L%VOGsFA>%tXt4mN}nA4QOQ)OV9Bk6n&-c~vpc*@^1fR9E*doC4a2Z)!vN*& zv-YppGV(;^nmLv6{SF9>JSW&hTZke+NZy5Qls7U z>O^{D*>LW81(!$5NZcFKs27jF^PL@6+laoZg}87(qWl-6V@d@|4-ay9D(|emx*w@z zbA4(6WYWBQ1}$SEZA}>QiE$)NXDX_&Hhq#kMEY^xy^q?Y_sQlqsNsir;fHjG0G%Fj@KdzNryc@G)$ltLY_KI= zKJrT^=E!<-Hwas&m9J4zh6$0++yT#M87z?0&XZ^H>!51k6-fGg7Kk7Pa2+;q#2q-7 z40FH7cXK!G+{u*a>sF@obWypO_ap4q#_IXqog56Ec zDd(AYY?gg|qKtqpZ^2Ti349bjX31I|$QhAr%kYM_3E>FjWk2@A_VjOHnu1iy>TGC8D!g|C^77YHO27QgbJm#eFi**xv)+)Y!#y>_-Y zC1yM}^Z)8T+~lL9fgtWJz&4*Vj!YS<{Pb{XOa2ZBo2xPwngSpH^d6|rVc3J}m9#Vo zhl!G+Vkjj>P*6~Pef`410wW_Mc+*IU`Pkvcm}DUFfoiRRw)Q0WMFOqy#*(f&!X`ArXAzQ`t8R1a_k%&&Mq#wT9qZ`!4oXu^@i=bqpEr0bTsnOO{nsfxAr z&SbS~b1;SH2y?Mq($wVS^|8{txtufznjTf=$efz|QMs@CTLVNe+NQG_K|xNc ze_D;YknWgS6@r}4;T7Umih}di#gWc_FPY;wi|OD~VFW8HYi=@vNvGy~%(8al>o?#v zg%t)JzhX})wiwKMmR~ghdciUvSO!yni4$7Z;MB&Di{q!m5OH!Udqq>Z{MgaaG4fSv z-%4kkfa_zO)!)tJYzIw>X|jx(W%;&AHy~Y0CG*iE2g0gK2Au(iOW}Z-nK(;XUunSo z!iGozUO44>bs_Fg(6SJ593L4Mm)o~;OL%#r@+jOb<8_imWJ&zV!H^Ce@CnW0L5GRj zBYej35LLK{xbHnaMMUAKSQYGB57kW$J5kv=z+W{+=>sO7Ckqkfwo91a2r@DO0hgle?CkP# zb;PIn`FVT>l_z+Zdv2##T;o8m@UWtS7chPq z{Y|BQUp9JHu-0Vz87JSU*r=Ejf3bITa}b->5^x0pyDzQdLbUsCk%qaJRxHrs{DlD7aX2nSzc(88PK}Ot<9x; zd&zEc5Pzr0kt>;2( zxIZDYcBR8c5E&meHFf*;O6~b>TYQyCmB;ej`%QFOVAB~&QCCSQny7&iNkeP+n5aMi zg0a2TMavFpM%73fZ6@rlz^^|!xVu|sW*Mip8xR_tubW46@m=$Gwyb5s~YTh^Cc z2Wv2AzEOz=)H=ml`h%)8`Wic43tB48Pf3W|DUP;zelRN5NKv4R zXgi5R&Gwc}Q=|`PjcoU*G2B?lo0}rxwOhW4!=G<9Qb6HxESxo9V<-ecv5&JdcNS#Q zJh-*>qquqg#SxJo{yBIppxM5^EC@|Y4cpuSy;-0kf?IEA8~tvg-95e_2@s8yRaIAi zKuT)1=>}XX=GW;D)zZc6SB!ltSWSK(U-1IsddkxGwq0q2rv*C1@UHDM(4oiqPBZB8 z1_lO3JpF`)h2>3Y^#+JAl$-kUVz6W!P@s8uxk^-y4pbEO41!z7QEjC?S zTSFE8Dd=7i8XC&e2mD;CU*iH)w3CxR-oU;o`N@^(Cy@~KxFQHaUJQD3pSOD?DK1vC z@)Q&l$ldpQKnM$>5E6XyB)zD}Zf~)JgwqUhD^aJ`Bd=;#s(EVnbS=+NU;nnZ)!9a= zBFLC!k76gc?M+5N{09x!)VDgGrKPk@G66^UEJV;W^Y0yzC!sAWxtcustzrB7`v+@7 zc{;TXR#PEen-Js~40Z21A}Y9sRk_ z@`z8v!^2@=Va&DXHe)5`&d$!@OlD?gtgWpPTv^OtW;LF${#uU;3yck`Ao z;?X@nN{q#;Ir;zT??>>tK@15!MwywJc$lB# zK39Y0>R5O)6I53Lmko6(fTL#~jynT{$)H25HFfxlRGM3G9%=7WA2)-%WRVpkp=sZ$ z(5PaWuB(NpX{9kYpd;%Wk3CRG88M)ek1)jT>?66WlNNR{_h;`Bnvb21lWo`hKXQ2E z*Ay1FVP8hn(v{<`T1qRs-=U53T58!#_nX|ZHz2cnFl532ZOW^ywEKL!CNlRs!3>Hk zjE(DZ@drz*3Re18_ig%e7jvafs$S!Q`r0P$7BtL%`3oJYlO5nN1^D}{<&wHs#WIHb zUozl8FReubEV)w3Mj&Vr-a*|a5f`qztLtXu6en)@i%)zmXUso|C;;jO$=-_xRz*IH zLG4}swvosYI1Y>1Ih8HVM!RBT?ID1zJx7w&>wiK0eRZ}XW-_AQ_h_)$zFuYrw5;#5 z;6u;FVB~D3qoO>wcg1(I*J9GwSQ>3}8OQ>pSx~?ZgRb#s-H!oI)5+tD)WOJ^k%KgRUE2 zW`UEsY97h!ThVgi@#bt;&3|2(Yf=m;JmI1MdWr?>HM)RtX}?e?_8-92G=Q@unCjhY zQnSQW1P>rP`y|!{6)WZOvUcM2}O|dF!S^RLXnFEs{4)<_no=M#aZPY0<$jV;@Byk#d_iTyx`yi$M;J_P@O)E zF!NP}^p%A!8g{%Z2>SG16%Ex|tPBW!i+vA4sxU?7b-<%GzTFUvbFZSWnb$KQ?pK3GoSQ7m@~?*S*2fK6E8JWA)Nlp@H{P18!O zY^k#JfEaIFt5FGpNZ_A8cLNg-U6)SpE`Ce(J?d{Tp6L+SS|O~P^9J=iOng3nC1bKN zHnul87Cv=2Hs=ir5K#AbNe5ZkeDjOOf^_#_iG3S%5jC=mu0YZFnVfTAND1_$X03%G zuMzkd5%r*BS58uZ9nSufMzn!{vNq>v8bMBY3;Kxzn_FR~d~x-5Z74%u>)kK5q0b$g zJ#8u|q33sDQ)(B>8K^C-oX#;X2bH|_?30L2=^!`jS4IBb&8R34HC%=1VuLPY$JvgU zwt1m3YewGH>>HOa&bh=ylZy zGL`uFi@;v^6W;8yo{w}VfAP%l%y+ib7C@vOS#AvsP07GIKZi>)FImNjbPwyfmIcYy z{n}VMTtLL;y>XeGV17QX3$QC3*}P;E21o=a-yE<0h$%GRXz~mFxDGM{VQ#R^;)~mm zR|-tM20%0IbX|^ikUZy^4Cu=U zRV5u^u}o5@-hietVHYeT2BcoB)`_KMoB7~!K&XY(iI7D=Bp>TIXqEs+g;agv_vni{ft5sEiXEGa!Fc)h_M|3@A~)-=(}7ViVk`P0~mL4b3(8psms?& z4h|Vs$exQPZlLVd%H=OrPJOP7;-;65v_7iMfBAX%HNn^puq3Z6*xsfATcs}GD5~Ny zapEdQL5jEA<+n4LiwUn}#R%Klc_hF{X`cU16B`76`3Y}QOpVae^^3ylwZnvuGGZ`w zfO#^qRM%D8^sL-nOdytQG753^uN#r^$yE#mD)EB%^wOQM>79-wd)cal}2pXSP=VyrYLD@gUtccq^`|1L+=z z&DqK;3uq1HFRz}LdQOcPsu3A>+=R^FhJtP_VBT8$cMUz-Gr@2fE$}`71|uU#>6Wsc zAulr6_SgGEdf9^A<0Ca32H7K7DYs+axSm$>zySwB_OaJHBL@#7IO9O~Rqbj_<`-zu z#*YsyojzfPm8xq1oXm@ih3Hebv3}Qu-LHz*<}m=iYLk;(j4!|V{0gpYp?28*^9j>T zqR4d|%?!13j9FPGkCOOYeO!@V2I0aRF3AJ*9$;gI9OpBNUFIIuby$OpkB)}sM{*Z} ztl>vEEkF~Dc95(s&ypvsp)v%|OzW zAxuUwJ-+u=L!`UbVq;}46OW4sf_}g!m`?^`8$gMmU)8_bMe*!-Max~nwvG~*5KNy9 z`+_cd{J2$ia^gT47$eaky?=3j@LGxo{Vwo5V_0niS@OjlBUWx&1}Vm@G^1;MFWr-> zHN{?`LS74SHpu=tse5eaxzk3*8N6eyTw`n6W$Z`#3>eP>tc<}`wG@-Xqe%3xOvC(( z6JsvU*j3_+Yc;Oe;45Gq<%)Z#1+L7z&1Lx@ZonnSS|{)dP-<2^Jt6wD?bZ zwt^RiI7lJ@H~pA17zLWngterSRHDwd)tn7zG_OWU&fk$Bk{V{8i@4k4S}3z><~&C;|EfHB~!GdzQ@|U`XXDXhGZP;UE=a2`i)!%<~*? zj|SWn$AWd{Y@jrsHTN!tl-q(F_%&=??a5kcI%N;9I?Q=YHcsXiBPPZyb)PkZ*g*-; zf7X!#?zHp#aJk1O+&Tl8FK}p3k7BcCAmuOway%bWFg?h#>`*UalgkovoODx2Zzeq2 zUt1WGfin>tCznOg{Yl7AK~A?%vR=vlR?Q-4vH8y`GHh=V1UPs?q5CoD=;oR69@NK) zOjr`4Ns1pH&&GMY)@GUoKpvuj*I63VQ4-o43s{e$nd?=jI%~#a0WFj;sWRXSUZ}9U z@+$gE9L+X`g06QjNv9SODdAc*Xp1G1%Oc-W3s1S(MS&s-$2~ekXUI&`!O_O)Xiyul zm9|;{LbSonH#={Fz|8|48-lHNE-fO^b1U4B^AO$)LcYQiX4T77Y8Qt(G;NFWYyOt) zTiG~W9ZrJ5rF&fX5LEUQXl;6x!_=;f-EK zO`zJd4aH1J^@Ehrz;|l(7Km>Oi?5jhSSx7~Irp-(Zo9YNXwd)wHi6%V&m-TzJzXkU zXc`qDjYI2$v$ZV17?J~W8M13zr(d+TZ)x7q?USomCDSf&d<2Up z17Vjv{~19n+(+VO&MR?lib?_#h=J08(tpGCw}b2f(?TlwV<ZUUL_a;++FD&*1x!Pq z4A;Wal0Oj+WLDR*no32n?mV*;+QhlMxt#CT>a~@IlZ*x6aaSYrC+P^Lz^`AU*o}K@ z+|L3+zkdC?GuLdWr?-6I5-tT(#mULZfVS|Bz7J)G)l<4vl|I|4o&eKyfuTLGwIjQN zTy%B6!J)%zS&Zq^+0jKjqQr)mib8DUZjS$Q=TMm^M)mRqBvqeI%R zGX2SrZaCGxE0P42maN?pwhEKJ;$iU^OL|E~F8fXrM2Xq^ClCH`O}Qx|KBcFp zbC?c3dh`gxTlX2!V4e;$#-DM(&vwOwza$8Mln5rT1|;9=s)>${j`hu(Uij(}>O5y> z?G*-r)+5Yy8_^lJtIsR^q5=LGPuKx3o5X_UxIRpU@#kO2!$88jf&C`M-7{380HYco zd;y)Zs%in%G%g}q&~+26DJ=ajlX&UhD@1^M9{KtCF)=ajXNP|q8!^1WDU?`D&|~~* z^t&%(X=-fDGWnF4l=P{J5|r=eT->2p`w}#5n&8-XY^N~gdWP=Hvd^y z-WCa>r$K!jVDY*?ewr9MPBWjGYmR?aRkHTsxW5-!FKh|v#pap>)}i_RpWFiKty^YOH4{@n)i0%_ zG>VKuc7{+W#~e)?GEy9F;e(UsFZ2-Mle8b36YAZ)|Q=l_pIn zDjM3?SaC#qFK%P4$T6+Pn>Rq+X)@4r$u(Kds)EAo6fH~5NyOX*2(Bm=LCAPA_>a0m zMO^%Eh5cG@U!P3qO>}J7AORQI&@h)8ZPD7*MaV)9srj&@xy;14I zr81OiPFJ)12^e<|4|8*J+&7Ac#AGm`U&cc#I-^b|>xSrHnCP}QWe5h&<74uZF?*19 z8-u`rV0QGsE=0}AVzf7ZW#E(!jXgB@`&}vp5Y66{KMx6&GnC$x zTHxP;dk(3k#YIJ+ONHT$03ks{3VQYsw48P2>bHXrqeiDM{Dd<+D}y)zGh5*eG!9y{ zOtiF1Tb@_0W@ct?Zs##Pwoy@0$qFmxr4t5jA7#^`xh&(PQogq>{QmuWHWk&2$-d&{ z-gH!wP<%Eg(*)@Ps52c+L=wsX9p{C0zZ_As{pXI3+m5)i^U3y1gG|E%HUlEAV0F+k znhO@vR#m0N_yc#cUZZ(o!Ch~lBjrGoIhz;0&6&oxaRp-j%AkU>Bl`bmH0b^dN)v#m z0^hj&8%6`&+Pxvk4QCr7WGv?ZnE})LncZ@CRaMnP$)Neg#VxK>Qc*m%sUnF8E?USAzNSUTpooJP%lGn6Ui9wT8CpQ#4gm~L^6&qDgo2+_fD;ke z-{?90;V@?dNX(%ltZhfz^U$eZ{~LL)&U&MaMkA->*Dn#K2=rUGd}C4T7FO&kc6V(e zUl9-);dX6DCqDG`1c-ZGK3e6f9Y-CE60tIY3n5&zje-%R12c8tMx+;$A>HTv-odw`F5vU=f z`ZeZu!5~}wRNoNHO2auo&`$?rsXemLY=l)9frkkPl!p%=g1a4EF+7TDYH$mxpr8P} zRqp3As)r97CSq+PQPh4L1a7t*G}rJdm}v6y{F|)I8Cp+FM2dF!3IrF63EuZ4F}(Jx zzr_8STqZS=dY*!=eP;~M-@ku@1rCQZmBAgj+tAP|7y_u)Ndto)j?qTh!E2>Qw1;3| z33*_))bXOnmuFjdyq`I&4==PuTzMQZeiqpoejof`s?eY#s#~GNvn{};} zF}DJ8*}9ficj5)qO_F@6;X810~#N&nyPX z)^)zJ9rQ>jF+dk$m+odw%=nb}j*sUGaMYiHKT87u3JjEZh#>+(?$x}xO?j1So;K*L zW#r^=sIo)Lw<5zcT~(3|hS(oO03N=^osy-qSsB(XMDtFY%Ffz4L@LFba(-#)>V)qO zCMLjEwsVK{kji_In7eU~uDu$4F(3iJK**okfGVs-5A>(JDYHQmmYXbj-y8&*(ZUEY zlJGuN53<6F>K`?pN|w=yo&42v+q+wa?#jUA54e}aZ5#ar82U46AyIDkQ#)aV;MMuw z(WI+f-q6)0)`JHR0GcAWMa9DM+psHcu@fEgd~J|oGJWO>uFY^UouQVcoNrWst2Ia$ zwt5Ns%WU(ST}4ll4o^ zf(*69hRHoYF~qICE9>+DfPqa?6js*OOc#UE%(_6oUe27`Vq8c9@d?x?Q?=ElNmXoD zD>6wvO3KP~c$h*$wX)|0P=&8bR!(~Qx3-1aDH$>^)ZP_4F#Ew8`Se*g^(D*DT{7W$ zTMT{H5r9Y8jC+$Cd~lZL=c_>Is`U(xh#30*SlX-->qfoEBz%=(4kFueDv^jSzoc0I8u#NqYn`ih$6G{UojiD@8Re_P zj%Y5&cSKN`31Hi0blmq}nO0gmgR7+9U%G3BpWg?3aXzcTMOVc4@87{3&Ck92C<(Rz z!g6`GIL}3BR@YNvcq=ezIoY0%hv2?j0tQ(|Y{2*l#1R~U{ck)%xDfB_&za_}L&jhj zEuih0+|m~L3sRzSEu_7xb^uxW)4`xeAnOwCjpv3VTJl^cYeiDZ(bhaHla7DgZ}P{I zHNqLfiqUBn6=DKS34)DwSN%S^oJp`WAW`0njPE5SJ-YWh6!@PeJc@T73F!4cRf>xl zr*La96jr`?*$8@&$Y6uKY;r9606!BH17P$gunl)6Zc6YQn6iC>R1Z95=Oul%`eJ2n zEvOi-Y`(qtlqMM0=hFphMkX(DEL?BG<7=7%sE34P&k>y#hCct{1Tb zHD%fJmB>M`2zEpa3lo?sJo7lM16L!FCHj@&R_D#k%)P7N+O0p;>%PAzpKR(%QCp-( z8iz_ptr?~2mxsSm-!qL$ZAsmTc0U_*tKCR4U;K_A{zfA*HJss>_nopCrYtY&wFys8 zQlGbY!g!lAlV@({_MHOFZT_=X;1YS~(#KO!9(#F?*{%OcWoS|RSnE;!hg zpa@6I=K^3OUcSB!R9hB&Crib~__osLL)B#Q)s@wwfKY^1;nBC;xHr`N-F!1Mn~p`s z*$NNrhX2F5A7>;h9tj0DnYvaw&V$yOH@lo#SMGBRMK5~5ycz{B+r z3jP+Q#eTS#{GKqkow29ZJT#ADlbk?T2^Plde!3zCm|@@sE&TTO_OGfM#U=-yS4W=j zn*z!mmVoi`?=;>Arz4MZPk}X!cAc41U0f_FMCicstjRrgUn3AHa*Xij+r=;@r4s$3 zR3r=nF|}c53@C%fB`MtJ&-o$#(1A^-#j_&xA6v?!3AL}U7DOkmSPIx+1mT2jNcabQ z30ZBMHOZ-8Fc?fcDtLI{=GSMnMX!sz2 zsQ-TVSvyZoTKWwrn?e;)wn^h@W7n=O9a6ANIHFX)NaRJ1`op-7KeX`9TMdmU-~HvD zxltJqalug$gs5ser+~_~r5-91@b9=hZvCiOD@9j2D*kA=*>oa6oB8h7-tBf!NNg*0 zGJW#!@SuLj8ksqMI_uBOY0|&^=kuMs+RF-CTNEqX5%JW@tMa{%8)@MKOdU6m&NHs6 z$)Yph-!<`EAAFdM*A5Xv#^<0^Ui>nY28>n|_`;JEtsN*}Z6 z8NS5>)3HX?2{;{Q+B#4vSrA-+B?q(B6&J7pw%8cViT?WzK2HGfW1_LWNKPN^m`xp( z#8#`(Q1j+2cHqX=M2$}X6Ni1(>oi6SBkfgWZ@I0ZPXi06BmmNZp#`Y(vGvCtadCTA z>=;zCF)%Q=JQ0@<)W zr`kNQ4Qgy`1UPbSbh}-!1iv>R4P~H~UzRARtHuzuAGKP$LHv zFfP@C{4tTTD1^F&BVcIHj#o=e?9-^aSWGW&jQeMKvDvw~U~&OoVCV@632PgD;DEuP z7#2CdTEn&{WWo##i~Q=yRh%S-i17Eu4&OKuI`lId%Du`n>QC=Mr*PY@N2k!!(@RcC zadvgpsg$#`+y8<^5=1UgVKR^jMrr_!#4SV%sTRG{X-b+}XXY3a6H6PYancSh`_=h}+)D|4{o zu5a_Zg`Qd?a507I_1>1>Y@#r*4r0)ZLp31DEE} zZf7roQ!ml%;o7Z?0Hma$y|b2cYw1HxI0(y4jvVHMfQS`G3^IE<8TdKjUQaHYTsFU6 z$%uGk`q1TGuvY4!NrnS}XHWznW2afo9Gdrp$erKYRumA=8Xg~pkUJKe{(>VgGz7XB z)5N(^j)I*0Rvks!Io|I5aLSzSjnTj;wU{{lOQ_R-s1wWQm@D`PvZ9i-@t1U~;1kFb zQ?ITyJ)XU-<&MRBEh}#t!5=W-2nI5F_84wh-6RRp{+VxQ`L#sXcZyUH?t0+WL=EbqUt)lx%4`o(;gMKlK&=4oLAzLZ-W zeJ=KvxXDO~9cTMR(tC{dU2~;McbpU>3~l6r_@F%Hj2(B~Sj{rbRpba4J3UnHbO^ zdG0PwvzEMYHkh*TBx@(6J}Zso6X;kyuLdm|`4PYK@3`@sf-#EO*@6hd-^U+vC+|`L zt;k;iQ|T|-P^$m=_oee=1Kb{?`M>o&o)mZ~J?Uh3faS7*9$q>P)nnUMRXTlwXkVg( zQMt`v_jK1g#iGDcoZ!n6_}rwsyK{D*r2wejt$D(K=I7mFdA!jsyoFk!y^ z3LukCTkf>u8oB#&ixKo$g=Zxmj;{i=+lW}f?pHxG#OgcainHgrTi0J!f8XVzD9?cP z{rOn)0c_Z$df)|S_L6PcYI$^FW1ULB+D)`Lvd2{B!vk7P11d__!S<;--|JDmz3bcx%<5OMHj<^L+ja{s(bqOn`Ajt56F_)C2y za@D;;Un35r*h-#gctC0+D|p%A){p8q?jf26zx{kmwJY`BjbDIF4e!c}9vn=6iepLK zxu&JL`Xe?bXPe5#YwE>}@Xp zsWhc<=o79)8V^FUnBP2BjMM~SbqO) z5Fa=*R!J;9(aFlDzdivj%W_e?tndGLPXacOP$_L!?<0!Dg`=X&-LEE@$oSb0e$$NC znnuxH`*=9IxdS*{oc+tU5^h7j&yjAI&T7K&0tDj|C=wj6;)kk;r`L?n)`GUQL_1(n=?A zBj7w2eBMfwe4i>WJ9`}vB_KN91_UCwtb~D4R({m-?;p3@u@xwTLFP3x%a&;1H?yLM zXnp!nXsr4`??VA^xvRe%ZKKOd(e+XSJHB1cpcW_NFbRu~*D*6&2kiyd{bgJtB5uHX zfyk{<-L=w}=7UW(TA<%1{1GjXn1h~)DHqVE0OX{n4SR*^Ycp+K)dzj$JYJvJh2c6e2}OIwQ#j$q_OM(DoO9-Sf*@JO z#l`#y+vok}csVG214+1KU%YriN5|u^t_H3-sOq@m2L_D>W7Vj}y?ON*lR4re3 zu{{#}_AJ2O0rm!r6XRfF7L0A=RUc{ycrrY9ABl-O^7+hX>n}q+YuF6{f$clZPspS* zeb+)81HR5dEs>JkbFji%UM)cgprS5nua2(hYbO&$eXQr2{Qd!= zA5eR%kJzSP>ChOCpS_-4u*vfMbo`*CeZ|$X*}w#DLwbdlIaBew9q0QN8TUz~^=TE; zfi{=(Z_qAkIvOdnPV&9`fVHu$Z7oygp+o>8WOsL$q`Vh{{RxGc&_y?4i$S!OmUwyv!yIGaIh>w)|-hG1n3qOyZj&o$NgwkyGM3&dq#?oIu^t`${{%KmQ z;R%M@0mzi6pV(6m1W@wHog9tg3I|hNU0qQ7({G`%TVea~S92buU`d7#Xom!P#P)Gf z|9gWBEmnd!Q9rf+pYFap9LxUyTXo->LWPoXyV?>Zku9r|nU$=Rm4qaFbw?DEGLl(l zW|8dCK$&GETiH8guixvuuI}&edwh@Qd47L9&tH$DuH$PuDok&Yh#{$ zp>3QgC@$W4>8R-wgOcMHMjs@&2~9V?-CQvgk-9d-_qyfb+(7N;7Tj$IeM%oiMUD6N z#`!+_Gx|1+@A>(8;3Lv!Q)qSFeKMNOE2`Eui5==o6G2P`hW`V$LTh!bcruu>qvM1= z!=WQb)EU#VGyVPjy$|loWDeuMo-}bS2qC5%rFP!1>iv&DeYZNfdhY#v1=oz$u+pky z{247)2G*tGJK94d!IRcLd2g+{xfp#u8G_tcGTqO;yqLQG@!5p>8&zw$-x_Cz%I{Cr z7#1u~F=2QIGYzvHtO`9Y4E`x{o)5Du3VrzSrNdMo;7%G->6|OE2_iS`kNA|5NGJ3R zRAu9H*?JcYGL384x>T;-k&XE>ug|P#`=U+A;Ot~tre;gmh**FnjeniB(YMAY)yDbl z7DsoUG+q;ODeEGLLNC4l$)7rZR~uttgd%EgKG3S=BtktE@2>cZ^55*M${))c#%O9q zrmE~ZT{NqUC#R-V8PjI%*7i;(N2(Qk6`gHUz4C5%?q#iPq(Z(B`H-5`I;;2*QvE}v z8Mcpk)~i$Jq@L@uu*6yHZ9_)w8Xq$a^}2Gdu5BqA$zMzxhNsDE(rhjpqmoeXP;+cJ z*FXB#V1560FLuTNB$+Ryel?phZ+ZP9;}sCoQ57l){c}j`5 z1ZiGoq>s_d|MAS?ilt>5S(fqbuyCz3ogHhJF{}(eqHeLgJliP~i+HbW%&UdjU1H8T z-Vk&=AX?4qa%axa1$_MY5wFpv-O})^Q?aLWxt3LK_b-8=CWUbsA|5GW1X<-pvwJ`3 z?AfzO;^mpuBN7^DAYSgfOLczQ!-L48aRbZ+IMk0)ghLM5`PD>;7f;sQyLYd`oex|7 zBEfP^j^H)bOfT2h!`8X#(vnMCifaa)l+;uy+STC8NtUfy#~b#9Re#O58y9t&@&1Ul zBUBV({qX3Qf7L4`You-qz)0Kk?KC69#E`fFs^kmx*xA_9$j*1O zIi#OZ1UHE!Im$&+9Nd8AlIa=v5(5z?uc%mLTaOqp7U&zoG8H^}>2xilGlx{;5nmWd z61duw*6uI9@(m(hyXq{5w5}FRknr)XEv#GV2kCxf-_VHCcq?Y+B^iMzT$-okh8#v( zE{`#hL|>ASkR$LhMHDLV>fW@cngEY*LkcX- zP2l^A3JRcGA75WYnHm7&xfXt#&gvNj6}9`m#1N2k;j#!Ne1u}DDSCRJA%?DG({3l_UpK8=kgYzk~N8Mj#N zd~c~LK$w;~=%-^P%2dfU>l*mPe&%_`!sNq&V%i7iDOVBXv^s_u)~U^QY!%$w=>23n zcCk)L7xPn;?m1K$ru8@*%!}fB^%1h={bi23-ccrrY&xgi$JY7Q{~70YA!o|a^@X?n z;Olb-Y+e51mL!|rRt7sHOHPDy9$bDeZWqaUkb(Y=;$m;le{*#XItg^GebB?$nQ3Jr zYkW=tO-E2bwd%E->XD~0j` zAhbVATo+eA3e&u@SA%}z#Bp=`6uWH3Ql&v4d^YPYJRr2?stUv3oeM%Bv`+j#}80dWkn*>4Lio znGSw(?U0oCb7N4;rT&-h`o_~6;|R;UWbmMed*;4aL$SGZ{+-(ye^!o&RWDt=gHbE5 zxl=)A0DZ1rUWFb(&i9LWt0`6s+G4rjluetinJ^8uI< z?V2hcesi_1P;U3gw<@=(k)kRA1d!dS;hi_cMEUFOoSkROtmY=MN-jNdLCq=^%_dK} zM4qq*x;qO`zGyRJ3V=l0Z8E*lWCublRD65!Og(c^t)q^VZHK!L-O$SU?KF8sO5CvV zx5l%#&(FFRB!4j9PTXExJ=mJ}DxFN(?#b9#p*P(&Z(X9sG+Hi^>7t(4h|dzZk{4(* zR`=51jwvjKw2i)}5F7DE0jO$t{Lx^Eh~_8Z&^Drm8wdrIsiZp~HlM z??xod%lUP3bOVXR;d+-$k<8I#l8Z(Igy^j-`O4O_;&y#ILZxE^{F{{gryp!2`%BT! zAg;1GMbiztyo58jLwwFW%xez#Ya-LGIw?|v^V`ZYf&F8~|ht_{3}g?Ib+Ck3SW)^%DH zL)YF5i5-6#@%t)bnr3}gAE6FiL&Kk?tc=pKuFQCC9g+ocQ(2(rV`6jWa6RLwNDNnP zZp<`pP6|Vikt%UJ_Jx|~`pJd9;O!NYJ!3~Yc5ja%>Zhrq(FgGMP(l85L!~&oYVhfD z50@DUO+tU{{pafNf_+nH!dI5_!Q-RfMuxrehqK7Up92A?oaGT%)!fNygG>$X=}l}?7EDiv-19SO67a;vs}korT6lFGwH4ZoBL5#dgxD2y1DFM1&=s< zY<0OZz0VwX>#hZ@ZJALai%(^JY&Apw=F0U)T7q;9w*~xVnNEZ27JcayQecDnqYv0z z+DFg|;RCP8E-U$$+WBozqdcG(LS{#3%6qZ1B2TrsNV0{ty@70n+J7dKHCUqE4oQ2V z$BkkOnX#h%-1oc{leH$4vrnlgP#n6G>pbb$s*UaB^#Te6MU@CbIhj`Ztu>5tyVbT$ z93mm_IWmhxe^#%^8%_Itk1+mlI*A-Uu|oGJktY2n_LaR8DH$T;HT<<(%oE5PR`Thf z$4}+)>%!CG+I2y_pW3wxL^);KMX==VQz{$2-xE_3n(Y~qs+b&j^0JiHrqQ%P0dkEh z5g}UtOztE7Up$r(tNU0o!=bM@mR5g$Jt3j*;h|iuN(}NB-ZYund7pt2OHP>i_4OtZ zB`cAdP}#}0_e8#Fk4gTyYJU2wO&=_medqq!Dyl~14c`ydFa^^t?tMjw1%+iAaS>*A z1f)Ji+(QL}$AqSQPVgff7bj-sMR*mfm2#dtR&Vd-!D{M%5`+ zw}+c_hAi;h7Y=$@-DC`NX?r9ptjOg3p9{^0t~)u$y`qKg{5OP4P(;x+cjM~K_N_L zmQ5UO3TC!n25KKwq}r0rRtv|1s8bZZ9L@2Z?pw-r`)47U|71To>Gyp6C zo&^F>r{f)~NyoypM}|nUWwHtN8&v^XTo zeQs{TFSx`-quaMHwx-`mgMl*eA4Na`G3sFY~*)7CVxr+HtI$F;%CY51uRBeDV(=@&Ju|M=@9p;kbgwX%*x%gl%J%d?S;RFZu(w3ENF?>V>j-$nJo3is$` zMJcb-PTA)vq&9=`(7B2Ri1+7%DUr5xo9-~BFl8jf+cayyC}h?am~AK=D}seL%lJ0Z9Wtg+dzhZyEoNnLu<+@6 zM~bNa3HW;ARjj*q-4uD<$Nl|iXw`)2G2X6IcyzU}wi0*S4OuO{qa$w?^Y12;N$Dpa zZtL!kVOlBRSI{a#(lB#f>+eevxYT88G2vn4l#1Z`FSw6W>+^V%h}315^A5?A+Q{#H zNYK2?KtuaWPc#cW^Zq`ecT;S`65>axT2 z7M^m9H2DixSVQWyKPY!5L83Vj%Aw}#XX$ANn$LA^+p zIomDLQFbyB{36qnBD%mbo8mRnM-eTeo9#Wek1oLj-jgb1CvbV za+<-2+MR$;n&&{i&mZ?0ylU=uwfj6~<^3@_F@9?OEoCyLxuY?#iazEbx~TRW=RjfIhiQdaNTls~cRT+PGGy~Qy9deX z852E5di^_Nev3jcSP2tT;)wEw8&lTz<#|?}<+WUISg9pnV*kBs44XJCI`qHB>2ktK zgCBWd+&tNqg%|}AXde`?*O!_*rYa`$uUz(ZB$o|RF}?S)Z;K7#XHhMCw5us@{_}(C zh6jX7iT--fBWD+~J-=032g};J5I0lgJ}0sYZuS;}nHz z4My2AkEAd{Ae1P?&^`Z(unLiL>s5((-#|k=yel3wdktJ3)V#lie*$^pXv#PR?n!!I zOwQ2{tfVX}aramD8g)DoX)xT$B+DT^y4Pe!B=fyKx`+GKw?@h2U$^UQjvav}1O7SbSvwr6aPXLXam zi+;Eb?}ntBf7@-hI~@P*f@G=3_C2}!grR2TvO}^h@$%zoIESR?Up{Uy{)?Ec4t#Zi zA+Ptj`;&F0^JM2701jh%OS1qjKEH&!Er%NeQBmZ%J z*>Hu(Y)exw+FWHzq^!yWLNce=Bi}oNL&6 z{>9(t={Fu7TFk!o{P(?gB)y*+>G2=hJSyo-Bjjp~Z&L}p^>S+9%6zG#Wy)6ai@el? zS~J#?S>^P&jZ`Ui&IMeMnc~*Y%gkn|CP1hRcOhvG>FNybo0Vkgc^9mpRmvcVY71O! z!-*6p7H1jWN&No29-oudHsk_khh2aRW7cJLM)8#v?DLH@U8ocQ8HnK#r}XYM6QuOg zKO{Sn9Y|@s`sX~`iA&qa!_Wn#+^LeJd_NN3e_L!4XX1^#I#o#_hPdwc>yUpvk;r#g z6v04T%9iAirdbi)Xh8LEpE-3oVu6PQPkKsvGcPS};>h+TK8_I|4ViRvC3iBSCJ}#W zL*Sq>SIziP<~M1$#dnGuhwi#BflE4DDNRdhKViY`M<$@wr19?!}ly8a)*}k>8Pzx)QAlg z*deLrb;=0;r4o?k_;4%e$#t3i$=>$ES$@t1zghHw%7z=dp<4iE!h7H^r}+pTAa0l? zZcrDI*jclVJRD)LY#w4N#McCfuW`r{i$wm~4a(~*cS!1JYN%j4DJzZ5%}JBf+>mQT zrgte^ZYtjV&vFWY5E|t*770p4JW*lkpAn5=9kNm?d%-|w5{pkIqqKx9wA#VQWVm6q%@%P9@v5*>9H5`90BqT+W z`W8+Cr;QE6_N#XXlOH~chyR(}9zbL$Qf2F_q9_wb2gjs`s)iZOtfabz&gdqy-;1aR zBJz~mQ(@&mUUzN0{MnB+@2q44iE*`IlhrR#$Zb2m`RCo8#NLan!Jp@5k9({k3YzgX z9Lb=bXyf^i@-9RKlS&T8;DY-B@YSZEO&q`J5OScAW$ULy{6xUw*}vVbnXh~qmJw;& zNOx~){{iBmR@%;O<&8U4fJdVGi6bO*Qi*X<_J@$pPpiLsTl9^aeD<86B9qhCUuNqi z?i^ow^BEHGzIWr8vM>sJsN5YtFi7fz=q0O=4{~yMB=EENp(S^k{6GBeJyoyxZFB}@ zQ!K)^RjsYU`}fBq<1H{N=IQCVek;pCE-v#g&t=H90m>B7`slO{2VbZ>(V$8~LIORi zkVz4`RO0^d;X}YIv#vFL9ddC+9Uap3g($3EUYun=c3E0SM+f;Gv8A6NFS^B6sfG31$W4VVDIhxJYO=$c_vGuzeOkT_TZxLGG7L$mWdnQ~vAB z%*>o56eS{c69KsBaxX!D?v6iJJbyN7>MNW(w>~s1G!uk60~!Z^#Q|=uSV`-lqrLscX#nO-xL7Z@K&O<IMnv+AlUWE#=-3&n`a$HS_F0|K7a52F8FLk;Oi;DfF{%?*G}n)8Rr^jNwrwyK9M zEid<1haDgxU8JR@)mZQeMHosELiY0OlxKQ+x?Xzr#Z zw*ixS@8eXhczXR{QnTRx{j#~#b2F`${yUHTNJ>TXw)+{%*nnjwRz6c@zqi&mZ(eIl zun4X0TUv6y^EO)K)h1-sFK~e>^z`&DriOJ+(dMjAZhRUQg$Mp}q+R~}XAUFQ8=o7^ z3JVGvYJTjEoXbCS?AVQ&k!Da5iTBn}11h;+FDx%FEv-z36ce=?=pc98xN5F*b-A;o zM1npSz+fbyGtC+~jReA~HI^5Ngm9TTN^%8$i8ioiHkPDQOQ zwOmVAcXLoi#?f(hEaK5463MX;SnL@_pWdENk5NMoE2i7JiD6VN3sJwqtl3+bR zOK;0Fba!{}nHCxKNYLSUeKFOtaQZ8d+7~ZgI4NPG66kbrl~U{-fEL1q$CeUKoNy_LZ8Dx5t#F*)fthUIo!ak%r= zs(!q0cz8EEyZ*7+Xzzmp59H#sXb9y|vMcfpUFt$>!tuP}v@oZ~ovw@3 z->qI>Z{_FZolH{DcU#W#ws&x#*?F{!5sgW~O2>2O6MF%9Uq$7BbkAcYvD`uagCw_w z`Y7tHEOY4gsxzHeZLD47oa>GDn~7y6WGxzf>c*7hT&>%Rsx&5Gx%BMqe<#T&-J>dg z-9fETXNJ+aQXf<~0fB`d60h6JQeG7eG0W-N>V$}a=FEO{Oa(|WFmzqINmH6G;BuP5 zS1w;Zuc|7u_weDvTwL-FucD*>f?PlrS4>PyNo?K_i-z{Hxu_n(bkwB*kzEAZt@E;) zk)D13>iE%6d2E4fOaR9S3O4w|GBY!2sHIAI4pPzlOQg{lZg8S%1pWVKUVCNBcXU3L z`EdCaaA%sCVU1`aZ6wQx}jI6GVucj4TE%V)FS-d=xdYVwuR0nSMOPXdj z_f?14?Qs#>yEpc-1FES&(4-&aQg{H)V@~@MFkg@wQDCnOWKx=L0>lmw-a<-*=}x8C z)tBv?M+Y<0YEO>Q4yz+=yqk~EmvkWHB=#29NmlVGrVW*E5u* z{y+H*tHk|jipDhz9*q!CR3$GjUmvf~sQ5>hu}+b*!|KX1L@s(2`}p`wO-!f=(|7gN z)zxifIeFCZUefx+a5NGO5Ted8g77*pf=G)QbGm=VM%>T_5`8#MpRZ-Lb@G( zqfs<>^zw%dVb$?p?bhUELu?ICFARz6eEt2<-GauV^GlUSVj>5;*j??>d09FZRg5TF zDOjFqotm74#q*vC@MM@GFtT)Vt>apDtx+=M4Mm~j0TPgSX7%yaS;D`-fRjU=t1Gi; zC_x{JI-;fdDS-L_zF&ooNx7Icv^ZW2Jv<9F-yh!N-cM9wcx>E`MaL9p$Wx3WMzg$> zY`e#EQ=s@&X<1pl>C+WxxZv$A8zX4{Gx}mPl@iRI_$6a8lzsXMqnL4iVd4JAmj>@S z0YQ5G`nA*S7=q|*WA_EpsA%eG#g~3Is4AhGsS@Ch&fa+)jdw3bm<3q93f||6)ovCwtD)-}CiNu}2bwx5XEUfi~>Ve)KweNt*ID^ap;5cxEr4n(0aR%~i*qU813 zOqNRbA(175F4_8_(MZ;X8JZK01_sZjc-PPcxGj5HqlODMwq;{dBTqE=vy5ua9;9E>D0A-C^`R##M3j@sro*l^waYh8_I@)TD&Z?sSiv6V4Z zl~*&HliJK)T?q87kovXOYVG4szhKl5)+%9hh;#;|r;o#z5Ji6t`gzzsL8qAl*J=GT zwnjo~D6>*_jsS_Bhhf+I^?uXx@^VDTM)xz&%KEDdO~cn+zI!r#P18M0f)y7bbgOQ7 zn)=RLbAC;?{}`1HK?31tpvu%Os%vT(5#>~!Z`Ptp8g9$Yz~Nnt)d!$_-@I6BAhfdNx`LAx=XQr*Krc#WL&#W&_a6AzCVQl-_Fm!B2{Halh# zJ9qd>gzh=PVCbcR6nCn$M6&!d{V4K%=^RFYLNq3NRELF!JCA-PWC9x9hh<(_!)6J0 z3|R7uFE5=6(uP|Ci7K-fcCdqtVDRin^UcvKThKJP!o3MJFx5!rokX809?#F+TUbsi z?^-#?!SOvw!Hbq%%x2)B_d);KwN@Fh-K(R-_2;{Pr~~H~!}X!8jB_&}%1i;zpFd|u z?QBv+@qFKIBZCj#^wU8MCATPZPOPb~E#7J$34ZeNj*99;BM(l!*R7)!wH6TmKKJjd z`L4H*FDMY<;o+grovIO5*st2F0C(}cvDDL}29l+ya!iGNN~(j`h+5(MO+0!9-=52C zrr#xGlQ>nkE|RD+tnf-inIflMB6+J5EgRW0!`*K26#ja&efR*gM@~)-sw7c>cC|4{ zo!PIfyu5$5)6H$B*#s?^c}PS{e}L*5Yc|-kptldu+EMT(VbkclOb4}F)YjH^q7M)d z_q9*xRaom$Q}@$ErwiVv z&FV0dgPuNh99;nl?)Iczbp^aRpa8pGt;m2z$kj9?c=fGPyPC7bdHG- zYlLR$=;$yqGD0hO>>$k6Wti0L4n6c)F>BaH(Z8P+_EV#odKurf&Z6bQ*^b4zZvSFv zRrQrz0#Of#3%iBl-t;Dv-n-Vto<#*;pOCn6Iz$g2Tza^}fBf->TCNmXq(~GdO+178WuN`IfQ;nBb5&n~^B$Wsg zvyH+{7GGVudHM2l?+aYPpk{q(YFaiH1_t$4N>{HA_D8uDp=RIL*SDmk1hE&9xvo2y zUj9((8_Ukf@bEV14z!4g61-`Q9yyS97^w?d+H=k{W1CT5s~8~^6cog;bLTf;G$8s3 zBlObIgTg3o=h7?4iH?@gQaE$wTh=vAL|706N=yO1)D~O>w@K4w+OlZgc(f(!fKue^ z3tWg%pwfC4<_+9o&T)9@j7Ki+CE*N)>8Ytj_#20vm(F1t8SS+ViFrmYS%k_)fww7x{#KE;4&j`>Lgj12nP@l!MXLN8AQJ zKf?@k^J~$pf79m8hl{76j*pNqOHuoummqn9+Vxzf)s2$dEux?$gZh6ULWAw&%twG@ ztc5jRDJso;w4*TJTV!!0GpN`Cr*U}G?@3$MR8diZ0hC%`qM{OzF$f_6nUFgZMYL%M znnf^!R&1^-D}UFPyO?tJ))iRO4P~fGbw51}{}#UQ)9~;D>s~n=*70sL&p{qeZL~L9 zm>Jd9(CF*#p8B}Q-)8vh5gs0b7%~4&m|Zwr#BQ{Ofa=cEoRgC)d%)`S{hcR-5b&$_ zdFRovKJCIQNofK$Z1BzlZjOF|{S64kWX<%Nj~{bT+5Py@Bk7*iY`X*w^;=FX<*R0o z6v_wohSD6U-p{uum8C1Tmr5&IfzIzT-Sq8-c21ubb=OXIx`Hy5&I^eWriT|g4k)6<{| zx&RzZ0hiW~Re0v~l%<*3ri~kcK41ty(+#+C`T6-k%R~Y5nUiypkB?7W+^y=#G5+E# z2p+g(@PUBZX*J{(rJ$bK#a-pm8Hd;Ce%=xx(19Sn|1tLLV}NiFNz>A5hH~jBa>+;R z4Pml<(dh2@sx#sYngmB>XmO|=4D-HdQHbX4$InWiQp-Xp@4c~=np&Xu6Zv?6Gm`u? zl5e4+xuR(R-d2OHk>_XhzNoIg4#xy7mRio=Sx?It@Ip1=VyedVeCtOZ-p&iZ&R$PR zOjKsZX&kYv=Eagv!oyR1MPNnS_TbyTvN6`Dvs83Bb_CGu<~J>8+OY$g$onIXsb~<- zv>WdrJ$m>sniavg=0Ybs`sM=+lguL?hl%+p)G6*Eeq3G&wuV1VoN?TpW~oQIVcB z_`qT+e;9U!;LZJE5o=3U*Q`a5`2jk@XiLlQSGBQu?UJ6MF0Za0RITlXb}2JnyWrE# z!Z5gX>()%1n9D+}&g-nK<;GobD!w*0-t&Hcb13A9HgE(2^bT`BpD8pII7~HMH&;XZ zum;5#I7SsRX!>*RoR4h}{9taqg2eDZ^aE|M_kvqpju?vchx+;{M5$&9v$C=(NRaIG;>0n%UcTy)V2_CvQb7)A2 z9eT*YNMdAS5;$>Ud0|VC!~i@G({VD-6RYb3?K1-0T33~FlX zB<)Jz_H1%S&{akz7Rw$VU-td%1T;il(Ryh&k_hp$gp2562Cv~9`Ns+G0|Ed+aHpA( z`NO4asEV1wDiBh5@+bz8rM{QhMMXt`g&`v#8W&i<1dE18HxDiGV%trBfq+;epKaS9 z^TCx}B`#yz285A?m6zaV!o!F1H1o;Mwv&{6hYm{6I7dhr6S5JkiI~wMc}5q|C%`)} zBt$>Y=2MX!&WJ(WkFWrpoYnbc-;wB2IHAwftK)YNFOeew18*UvCK>oimwdi=n? zeQ`-iX74;UZf4~A`SWLDQc?xlnbg+OZrG5RlCp%j0y>G9R6otK8`rxfM#?b#@*H2o zUSYX~Q$*~>+7RVHLV{h)`Q-8A_2&HGB18(?+H~kA@I68IqGd@*N#r8NItmxjN~x!} zS4LKrdCT3|@lF$CW6suKi*GM)Gi_KCU4LyIm4s5fVgx4%C)0f9cb%U!>iHNj6mmsQ z?tY89h6W>!sF@la92&}ll?4HJ&F2CaBmKs$TW#$-(XPX0aFS=0cDSd<@YMw*eHvmm zAc&CkAPUSZWIgetR4yjNEE}%C=VvOcBHU3gUJw{KU~*1k4n!l{c0qxga$tdZC^ zpl?8?7TrW_#yecVs5qu5A&L+Ew5-gA2Q(Vcg;N%Uwmdg+1+hUXDXFYcoPPuN!ELI~ z&&%sHdR@%K7D2`lp;ZK7;F|$j`R%0vq%{H<{6Uo!6-mj-2zFk&6#ZTs=s?80auH3* z&4m>F4uBdkhUfEuCW5(ASXhYYWBh|+*egiRI86`15E~dD=OOvCi?6_02cl@PYu2n8BE?o8=pGz=QC&u|h8}~f2nkgiB~FXO zn##*f5c@?0eWHIe)qolydQeWFYR*X{qN#*r{`!OX9YG>ztlbsq$+2aLgnz#Gha_?N v^C5y4&?tAuZF1Dp<~{@_e*XWueR(78@Coz7N<-zu+|J4 Parameter estimation using proteomics and flux data · DifferentiableMetabolism.jl
GitHub

Parameter estimation using proteomics and flux data

using DifferentiableMetabolism
 using AbstractFBCModels
-using Symbolics
+using FastDifferentiation
 using ConstraintTrees
 using COBREXA
 using Clarabel
@@ -19,7 +19,7 @@
 delete!(model.genes, "g3")
 model.reactions["r4"].gene_association_dnf = [["g2"]]
 model.reactions["r1"].lower_bound = -1000.0
-model.reactions["r2"].lower_bound = -1000.0
-1000.0

now models looks like this

simple_model

kcats = Symbolics.@variables r3 r4
+model.reactions["r2"].lower_bound = -1000.0
-1000.0

now models looks like this

simple_model

@variables r3 r4
 
 reaction_isozymes = Dict(
     "r3" => Dict(
@@ -40,9 +40,9 @@
 
 gene_product_molar_masses = Dict("g1" => 20.0, "g2" => 10.0)
 
-Symbolics.@variables capacitylimitation
+@variables capacitylimitation
 
-true_parameter_values = Dict(capacitylimitation => 50.0, r3 => 2.0, r4 => 3.0)
+true_parameter_values = Dict(:capacitylimitation => 50.0, :r3 => 2.0, :r4 => 3.0)
 
 km = build_kinetic_model(
     model;
@@ -85,12 +85,17 @@
         bound = nothing,
     )
 
-estimated_parameters = Dict(capacitylimitation => 50.0, r3 => 5.0, r4 => 1.0) # initial values
-parameters = [r3, r4] # will differentiate against these two parameters
+estimated_parameters = Dict(:capacitylimitation => 50.0, :r3 => 5.0, :r4 => 1.0) # initial values
 η = 0.1 # learning rate
 
 losses = Float64[]
 
+kmKKT, vids = differentiate_prepare_kkt(
+    km,
+    km.loss.value,
+    [:r3, :r4, :capacitylimitation],
+)
+
 for k = 1:150
 
     sol2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(
@@ -104,16 +109,14 @@
 
     push!(losses, sol2.loss)
 
-    sens, vids = differentiate(
-        km,
-        km.loss.value,
+    sens = differentiate_solution(
+        kmKKT,
         x_vals,
         eq_dual_vals,
         ineq_dual_vals,
         estimated_parameters,
-        parameters;
     )
-    measured_idxs = [1, 3, 11, 12]
+    measured_idxs = [1, 3, 12, 11]
 
     x = [
         sol2.fluxes.r1,
@@ -125,14 +128,14 @@
     dL_dx = x - measured # derivative of loss function with respect to optimization variables
     dL_dkcats = sens[measured_idxs, :]' * dL_dx # derivative of loss function with respect to parameters
 
-    estimated_parameters[r3] -= η * dL_dkcats[1]
-    estimated_parameters[r4] -= η * dL_dkcats[2]
+    estimated_parameters[:r3] -= η * dL_dkcats[1]
+    estimated_parameters[:r4] -= η * dL_dkcats[2]
 end
 
-lines(losses; axis = (xlabel = "Iterations", ylabel = "L2 loss"))
Example block output
estimated_parameters
Dict{Symbolics.Num, Float64} with 3 entries:
-  r3                 => 6.88387
-  capacitylimitation => 50.0
-  r4                 => 0.993201
true_parameter_values
Dict{Symbolics.Num, Float64} with 3 entries:
-  r3                 => 2.0
-  capacitylimitation => 50.0
-  r4                 => 3.0

This page was generated using Literate.jl.

+lines(losses; axis = (xlabel = "Iterations", ylabel = "L2 loss"))Example block output
estimated_parameters
Dict{Symbol, Float64} with 3 entries:
+  :r3                 => 2.05711
+  :capacitylimitation => 50.0
+  :r4                 => 3.00787
true_parameter_values
Dict{Symbol, Float64} with 3 entries:
+  :r3                 => 2.0
+  :capacitylimitation => 50.0
+  :r4                 => 3.0

This page was generated using Literate.jl.

diff --git a/dev/index.html b/dev/index.html index d3348ed..4a073a2 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -README · DifferentiableMetabolism.jl
GitHub
+README · DifferentiableMetabolism.jl
GitHub
diff --git a/dev/objects.inv b/dev/objects.inv index 12c98cf10af69e4216da6617b25fcc3497306d26..90742aab9e241e8e9edaa681353748703bc99e38 100644 GIT binary patch delta 1136 zcmV-$1dsd13G4}wQ2{fNQa*o8Q{y-gz57?FYBx6oJ3uW})LudYRZ~06?qt}xx|CQD z4Yrk$WC&3F_iahGB`5JGA2Sz7WcBOr*WFKYSQ^nak}5il_yXY%2+hvN)64VV16~n9 z5~>7S;Dkpsmj1RKj}Oqljh8X;IGa~{JZyE{0<5xalBrh=x55k)ag zA0HBGv3yTW`wr05d&bfmWS=_syF4Eiu~UJmR4f5gp87s_?N_d#U2n&Qox zwnVe5MIeY?O;Qqm#xs9lwSw4ID<+0wVqHwEi&@31 z=@OS?uRi(2No2f99^X++uw=<9_6MSkS&d*o4FZB_4T9w|7)OzdS7ODnTJ}K)d*tkq z;t1?>?7x9>RnrC)*7J=;c_*#`l%UYCYeiEDWHn(pBczRxGYWss*DFD#{*37zG(cWc z#bTxyZ31vT}m2(z2W;Pxe=jGo7tmQowl$m z+l$~U+Lx<+tG|DaIh1^N5FL@Fow(087U2^&y;UjhCc4B8Z&i3@@Xrbqiapv<0d<rv+Ca z1#(~i-CA_tpZCz_LsPOn^o=eQ7f)Xw67c-|jK0UGz}0_}N2rqa&_G8gXS+h(vFyH& zPNt*7$>{Kpogv#tHf(^gfoqeKX%nF2;mD5sxL8RIc-2KncC<7_p#0U23Jai>aR}JnWzeM*m#qq}dWAr+_wYRc- zmj6qE%qPpjc&;%Q>f&s2X%3V9%@s)8Zk>T@-Y>xmn*DV;qa^OlMQt#UZ9ztRKxOuY z$;_%EteEiN|7>(L86AyHkPYN;d&1GkFH3mO)G~kcu4Y5qWy%&RKX3NXMa+r%^LW6) z(<8_i^^y+Poik=B54!!(skAg_wI&vkblv!X^rrw&CIJFPW}0Li!-S?Vdo}H2)w1N!#i`ibrWtn zeqdzVH>2kxUv_toOT>H6?i)|;9T8?JrdW zTiZ${83+vj`>iG0k`w>t$MF(~wA$V8+uhG{!X$t@B2{n_aG%f*ge=ZyXBX#_4|L50 zi$RIdr!fy8lCD|J_5u03Tp~VCq!KJl6kRe20l{^)&YPv!mEEP5Rj;lDWH(t zZ(=aM+(>P92jtAz3zKVNc5V2&$VWx)ASg%`#(0z`Zp?a%5+DfToJ(aEcrCaB9*4dp zED6X}l-`kmDK-(H6ZF9iP%a*jMe*-s3rm62!9q(cZb>4xM9R9T`t98)vTRN1c1=67 zscR8TVo;NmgrD*hS*?Hgw;ZBZYE8c&5GKB}D209poDx%AOlpfsv6wU~Ce4aj_1DlP z7UQ5k`L>BQbV=b(K_c;pY{|bgC)RxOg|k3sEO~STlu;SR>&P9duIgkSbYj{ALB!gF zsmEj%1QxTTic_^3;vLTud!8sx@O?)dXAhC=kjvqG z#gc2#z#OORmNP$W6M<>F_rD5Z-^5OhH@7n5%o zC+p`ev~x&ys7o>17<=XZq5WcE!9%zM^uuu6M@kiJ|1X zhrEb1>BW6+Sfqc??2J`Kx1YQ+J7QIdm6N~Ls72wu5jxbNWdt!@rePGc<9E-I{2VFH zxpd>kMXum7-v6ou|fp2k7HUZ$uy5663CmTVr$kLU!^tav~wi4AZy{djsj zJ)BPu-=o9n$FsIXs*g~UiIs~>R13MEeeNvU?~i+AdsBb+%O3d-zT)EHQ^cA2{rIuZ z!dd96hf?kLzozjX*)gx|KMtp?1~hV32fJq3en)Bz&9(5&mdP01@>Ql+TwP5GjP-T5O~fYdAFfjd5&g}ZIZUq4zJ;WT%T0iUaYs%TI%WfuarfXJj1Fs zYalx++US4c8k7gw=}>K4z@v*XN^07(coNOsN@Xv{pl$(@5iIUJoW~UcZ)9(x6m5CD za541r(EDlbruoC?-vpnS4~B)=GV|%HPp9(>{jlEOUP0i-bb-TvyTT}B@yE#mSTvZ6 z#$*s}K}L5%CHmsX!b~=<==AXKVtO>69?f);-7U_iM|;D-MWwtAK$Am z30M^2=mfBq-Y=0?V$=*aILo@NznZ!4@2?<7_WSDzkJA)A<>grwb~whD5Y~3UHjEcD ze3iNLZuA|_wx=%F$eOi#&-=A3P0tz)yplXi2z!|E&u=k diff --git a/dev/reference/index.html b/dev/reference/index.html index ce685c7..3b07fc4 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -1,79 +1,75 @@ -Reference · DifferentiableMetabolism.jl
GitHub

Reference

Parameters

ParameterLinearValue

DifferentiableMetabolism.ParameterLinearValue — Type
struct ParameterLinearValue <: ConstraintTrees.Value

An extension of ConstraintTrees.LinearValue where the weights are parameters.

ParameterLinearValues can be combined additively and multiplied by real-number constants.

Multiplying two ParameterLinearValues yields a quadratic form (in a ParameterQuadraticValue).

Fields

  • idxs::Vector{Int64}

  • weights::Vector{Symbolics.Num}

source

ParameterQuadraticValue

DifferentiableMetabolism.ParameterQuadraticValue — Type
struct ParameterQuadraticValue <: ConstraintTrees.Value

An extension of ConstraintTrees.QuadraticValue where the weights are parameters.

Behaves similarly to ConstraintTrees.QuadraticValue. Thus, the cleanest way to construct a ParameterQuadraticValue is to multiply two ParameterLinearValues.

Fields

  • idxs::Vector{Tuple{Int64, Int64}}

  • weights::Vector{Symbolics.Num}

source

ParameterBound

DifferentiableMetabolism.ParameterBetween — Type
mutable struct ParameterBetween <: ConstraintTrees.Bound

Representation of an "interval" bound where the lower and upper bound values are parameters. Since Symbolics.Num is a subtype of Real, the bounds could also be any real number, but they are converted by the constructors to Symbolics.Nums.

Fields

  • lower::Symbolics.Num

  • upper::Symbolics.Num

source
DifferentiableMetabolism.ParameterEqualTo — Type
mutable struct ParameterEqualTo <: ConstraintTrees.Bound

Representation of an "equality" bound, where the bound value is a parameter. Since Symbolics.Num is a subtype of Real, the bound could also be any real number, but it is converted by the constructor to a Symbolics.Num.

Fields

  • equal_to::Symbolics.Num
source

ParameterIsozyme

DifferentiableMetabolism.ParameterIsozyme — Type
mutable struct ParameterIsozyme

A parameterized isozyme struct which includes parameters in the kcat_forward, and the kcat_backward. If the reaction does not have a turnover number,nothing can be used.

Fields

  • gene_product_stoichiometry::Dict{String, Float64}

  • kcat_forward::Union{Nothing, Symbolics.Num}

  • kcat_reverse::Union{Nothing, Symbolics.Num}

source

Parameterized models

Kinetic models

DifferentiableMetabolism.build_kinetic_model — Method
build_kinetic_model(
+Reference · DifferentiableMetabolism.jl
GitHub
Reference
Parameters
ParameterLinearValue
DifferentiableMetabolism.ParameterLinearValue — Type
struct ParameterLinearValue <: ConstraintTrees.Value
An extension of ConstraintTrees.LinearValue where the weights are parameters.
ParameterLinearValues can be combined additively and multiplied by real-number constants.
Multiplying two ParameterLinearValues yields a quadratic form (in a ParameterQuadraticValue).
Fields
  • idxs::Vector{Int64}
  • weights::Vector{FastDifferentiation.Node}
source
ParameterQuadraticValue
DifferentiableMetabolism.ParameterQuadraticValue — Type
struct ParameterQuadraticValue <: ConstraintTrees.Value
An extension of ConstraintTrees.QuadraticValue where the weights are parameters.
Behaves similarly to ConstraintTrees.QuadraticValue. Thus, the cleanest way to construct a ParameterQuadraticValue is to multiply two ParameterLinearValues.
Fields
  • idxs::Vector{Tuple{Int64, Int64}}
  • weights::Vector{FastDifferentiation.Node}
source
ParameterBound
DifferentiableMetabolism.ParameterBetween — Type
mutable struct ParameterBetween <: ConstraintTrees.Bound
Representation of an "interval" bound where the lower and upper bound values are parameters. Since Expression is a subtype of Real, the bounds could also be any real number, but they are converted by the constructors to Expressions. 
Fields
  • lower::FastDifferentiation.Node
  • upper::FastDifferentiation.Node
source
DifferentiableMetabolism.ParameterEqualTo — Type
mutable struct ParameterEqualTo <: ConstraintTrees.Bound
Representation of an "equality" bound, where the bound value is a parameter. Since Expression is a subtype of Real, the bound could also be any real number, but it is converted by the constructor to a Expression.
Fields
  • equal_to::FastDifferentiation.Node
source
ParameterIsozyme
DifferentiableMetabolism.ParameterIsozyme — Type
mutable struct ParameterIsozyme
A parameterized isozyme struct which includes parameters in the kcat_forward, and the kcat_backward. If the reaction does not have a turnover number,nothing can be used. 
Fields
  • gene_product_stoichiometry::Dict{String, Float64}
  • kcat_forward::Union{Nothing, FastDifferentiation.Node}
  • kcat_reverse::Union{Nothing, FastDifferentiation.Node}
source
Parameterized models
Kinetic models
DifferentiableMetabolism.build_kinetic_model — Method
build_kinetic_model(
     model::AbstractFBCModels.AbstractFBCModel;
     reaction_isozymes,
     gene_product_molar_masses,
     capacity
 )
-
Essentially an enzyme constrained metabolic model, but the kcat can be a arbitrary symbolic function.
source
Solving models
Solving models
DifferentiableMetabolism.optimization_model_with_parameters — Method
optimization_model_with_parameters(
     m::ConstraintTrees.Tree{ConstraintTrees.Constraint},
-    parameters::Dict{Symbolics.Num, Float64};
+    parameters::Dict{Symbol, Float64};
     objective,
     optimizer,
     sense
 )
-
Construct a JuMP model by substituting parameters into the model, m. Set the objective and the optimizer, as well as the sense similar to COBREXA.optimization_model.
Converts all inequality constraints to the form A * x ≤ b.
source
DifferentiableMetabolism.optimized_constraints_with_parameters — Method
optimized_constraints_with_parameters(
+
Construct a JuMP model by substituting parameters into the model, m. Set the objective and the optimizer, as well as the sense similar to COBREXA.optimization_model.
Converts all inequality constraints to the form A * x ≤ b.
source
DifferentiableMetabolism.optimized_constraints_with_parameters — Method
optimized_constraints_with_parameters(
     m::ConstraintTrees.Tree{ConstraintTrees.Constraint},
-    parameters::Dict{Symbolics.Num, Float64};
+    parameters::Dict{Symbol, Float64};
     modifications,
     objective,
     optimizer,
     sense
 )
-
Solve a model, m, by forwarding arguments to optimization_model_with_parameters. 
Optionally, set optimizer attributes with modifications. If the model does not solve successfully return nothing. Otherwise, return a tuple of the solution tree, and vectors containing the values of the primal variables, the equality constraint dual variables. 
These duals are ordered according to the constraint output of calling equality_constraints and inequality_constraints respectively.
source
Pruning models
DifferentiableMetabolism.prune_gene_product_molar_masses — Method
prune_gene_product_molar_masses(
+
Solve a model, m, by forwarding arguments to optimization_model_with_parameters. 
Optionally, set optimizer attributes with modifications. If the model does not solve successfully return nothing. Otherwise, return a tuple of the solution tree, and vectors containing the values of the primal variables, the equality constraint dual variables. 
These duals are ordered according to the constraint output of calling equality_constraints and inequality_constraints respectively.
source
Pruning models
DifferentiableMetabolism.prune_model — Method
prune_model(
     model,
     ec_solution,
     flux_zero_tol,
     gene_zero_tol
 ) -> AbstractFBCModels.CanonicalModel.Model
-
Prune away reactions, metabolites, and genes from a model using ec_solution, which is the result of an enzyme constrained kinetic simulation. Fluxes and gene product concentrations smaller than flux_zero_tol, gene_zero_tol are removed. Metabolites that do not take part in the remaining reactions are also removed.
source
DifferentiableMetabolism.prune_reaction_isozymes — Method
prune_reaction_isozymes(
+
Prune away reactions, metabolites, and genes from a model using ec_solution, which is the result of an enzyme constrained kinetic simulation. Fluxes and gene product concentrations smaller than flux_zero_tol, gene_zero_tol are removed. Metabolites that do not take part in the remaining reactions are also removed.
source
Differentiation
Differentiation
DifferentiableMetabolism.differentiate_prepare_kkt — Method
differentiate_prepare_kkt(
     m::ConstraintTrees.Tree{ConstraintTrees.Constraint},
     objective::ConstraintTrees.Value,
-    x_vals::Vector{Float64},
-    eq_dual_vals::Vector{Float64},
-    ineq_dual_vals::Vector{Float64},
-    parameter_values::Dict{Symbolics.Num, Float64},
-    parameters::Vector{Symbolics.Num};
-    scale
-) -> Tuple{Any, Any}
-
Differentiate a model m with respect to parameters which take on values parameter_values in the optimal solution with respect to the objective. The primal variables x_vals, and the dual variable values eq_dual_vals and ineq_dual_vals need to be supplied. 
Internally, primal variables with value abs(x) ≤ primal_zero_tol are removed from the computation, and their sensitivities are not calculated.  
source
Internals
Constraint tree extensions
Optimization problem builders
DifferentiableMetabolism.equality_constraints — Method
equality_constraints(
+    parameters::Vector{Symbol}
+) -> Tuple{Tuple{SparseArrays.SparseMatrixCSC{FastDifferentiation.Node, Int64}, SparseArrays.SparseMatrixCSC{FastDifferentiation.Node, Int64}, Any, Vector{FastDifferentiation.Node}, Vector{FastDifferentiation.Node}, Vector{Symbol}}, Any}
+
Prepare a model m with objective for differentiation with respect to parameters.
source
Internals
Constraint tree extensions
Optimization problem builders
DifferentiableMetabolism.equality_constraints — Method
equality_constraints(
     m::ConstraintTrees.Tree{ConstraintTrees.Constraint}
-) -> Vector{Tuple{Union{ConstraintTrees.LinearValue, ParameterLinearValue}, Symbolics.Num}}
-
Return all the equality constraints of m as a tuple ({Parameter}LinearValue, value) representing {P}LV == value for each entry.
source
DifferentiableMetabolism.get_equality_constraints — Method
get_equality_constraints(
+) -> Vector{Tuple{Union{ConstraintTrees.LinearValue, ParameterLinearValue}, FastDifferentiation.Node}}
+
Return all the equality constraints of m as a tuple ({Parameter}LinearValue, value) representing {P}LV == value for each entry.
source
DifferentiableMetabolism.inequality_constraints — Method
inequality_constraints(
     m::ConstraintTrees.Tree{ConstraintTrees.Constraint}
 ) -> Vector
-
Return all the inequality constraints of m as a tuple of bounds converted to the form ({Parameter}LinearValue, upper) representing {P}LV ≤ upper for each entry. 
source
Symbolics extensions

+

Return all the inequality constraints of m as a tuple of bounds converted to the form ({Parameter}LinearValue, upper) representing {P}LV ≤ upper for each entry.

source

Expression evaluation utilities

diff --git a/dev/search_index.js b/dev/search_index.js index 1b0cfce..71a1a7e 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"EditURL = \"3-parameter-estimation.jl\"","category":"page"},{"location":"3-parameter-estimation/#Parameter-estimation-using-proteomics-and-flux-data","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"","category":"section"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"using DifferentiableMetabolism\nusing AbstractFBCModels\nusing Symbolics\nusing ConstraintTrees\nusing COBREXA\nusing Clarabel\nusing Tulip\nusing JSONFBCModels\nusing CairoMakie\n\ninclude(\"../../test/simple_model.jl\") #hide","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"prune model for brevity","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"delete!(model.reactions, \"r5\")\ndelete!(model.genes, \"g4\")\ndelete!(model.genes, \"g5\")\ndelete!(model.genes, \"g3\")\nmodel.reactions[\"r4\"].gene_association_dnf = [[\"g2\"]]\nmodel.reactions[\"r1\"].lower_bound = -1000.0\nmodel.reactions[\"r2\"].lower_bound = -1000.0","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"now models looks like this","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"(Image: simple_model)","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"kcats = Symbolics.@variables r3 r4\n\nreaction_isozymes = Dict(\n \"r3\" => Dict(\n \"isozyme1\" => ParameterIsozyme(\n gene_product_stoichiometry = Dict(\"g1\" => 1.0), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = r3,\n kcat_reverse = nothing,\n ),\n ),\n \"r4\" => Dict(\n \"isozyme1\" => ParameterIsozyme(\n gene_product_stoichiometry = Dict(\"g2\" => 1.0), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = r4,\n kcat_reverse = nothing,\n ),\n ),\n)\n\ngene_product_molar_masses = Dict(\"g1\" => 20.0, \"g2\" => 10.0)\n\nSymbolics.@variables capacitylimitation\n\ntrue_parameter_values = Dict(capacitylimitation => 50.0, r3 => 2.0, r4 => 3.0)\n\nkm = build_kinetic_model(\n model;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nsol, _, _, _ = optimized_constraints_with_parameters(\n km,\n true_parameter_values;\n objective = km.objective.value,\n optimizer = Tulip.Optimizer,\n)\n\nsol.fluxes","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"create a loss function","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"measured = [\n sol.fluxes.r1,\n sol.fluxes.r3,\n sol.isozyme_forward_amounts.r3.isozyme1,\n sol.isozyme_forward_amounts.r4.isozyme1,\n]\n\nkm *=\n :loss^ConstraintTrees.Constraint(;\n value = 0.5 * (\n ConstraintTrees.squared(km.fluxes.r1.value - measured[1]) +\n ConstraintTrees.squared(km.fluxes.r3.value - measured[2]) +\n ConstraintTrees.squared(\n km.isozyme_forward_amounts.r3.isozyme1.value - measured[3],\n ) +\n ConstraintTrees.squared(\n km.isozyme_forward_amounts.r4.isozyme1.value - measured[4],\n )\n ),\n bound = nothing,\n )\n\nestimated_parameters = Dict(capacitylimitation => 50.0, r3 => 5.0, r4 => 1.0) # initial values\nparameters = [r3, r4] # will differentiate against these two parameters\nη = 0.1 # learning rate\n\nlosses = Float64[]\n\nfor k = 1:150\n\n sol2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(\n km,\n estimated_parameters;\n objective = km.loss.value,\n optimizer = Clarabel.Optimizer,\n sense = COBREXA.Minimal,\n modifications = [COBREXA.silence],\n )\n\n push!(losses, sol2.loss)\n\n sens, vids = differentiate(\n km,\n km.loss.value,\n x_vals,\n eq_dual_vals,\n ineq_dual_vals,\n estimated_parameters,\n parameters;\n )\n measured_idxs = [1, 3, 11, 12]\n\n x = [\n sol2.fluxes.r1,\n sol2.fluxes.r3,\n sol2.isozyme_forward_amounts.r3.isozyme1,\n sol2.isozyme_forward_amounts.r4.isozyme1,\n ]\n\n dL_dx = x - measured # derivative of loss function with respect to optimization variables\n dL_dkcats = sens[measured_idxs, :]' * dL_dx # derivative of loss function with respect to parameters\n\n estimated_parameters[r3] -= η * dL_dkcats[1]\n estimated_parameters[r4] -= η * dL_dkcats[2]\nend\n\nlines(losses; axis = (xlabel = \"Iterations\", ylabel = \"L2 loss\"))","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"estimated_parameters","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"true_parameter_values","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"This page was generated using Literate.jl.","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"EditURL = \"1-parametric-models.jl\"","category":"page"},{"location":"1-parametric-models/#Parametric-constraint-based-metabolic-models","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"using DifferentiableMetabolism\n\nusing Symbolics\nusing ConstraintTrees\nusing COBREXA\nusing Tulip\nusing Clarabel","category":"page"},{"location":"1-parametric-models/#Load-and-solve-a-simple-model","page":"Parametric constraint-based metabolic models","title":"Load and solve a simple model","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"The code used to construct the model is located in test/simple_model.jl, but it is not shown here for brevity. Below is a visualization of the model.","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"include(\"../../test/simple_model.jl\"); #hide\nnothing #hide","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"(Image: simple_model)","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Build a basic ConstraintTree model without parameters","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"m = COBREXA.flux_balance_constraints(model)","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Solve normally","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"base_model =\n COBREXA.optimized_values(m; optimizer = Tulip.Optimizer, objective = m.objective.value)\nbase_model.fluxes","category":"page"},{"location":"1-parametric-models/#Add-parameters-to-the-model","page":"Parametric constraint-based metabolic models","title":"Add parameters to the model","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Make bound of r2 and mass balance of m3 parameters","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Symbolics.@variables r2bound m3bound\n\nm.fluxes.r2 =\n ConstraintTrees.Constraint(m.fluxes.r2.value, -2 * ParameterBetween(r2bound, 0))\n\nm.flux_stoichiometry.m3 =\n ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2)","category":"page"},{"location":"1-parametric-models/#add-parametric-constraints","page":"Parametric constraint-based metabolic models","title":"add parametric constraints","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Symbolics.@variables p[1:4]\n\nm *=\n :linparam^ConstraintTrees.Constraint(\n value = p[1] * m.fluxes.r1.value + p[2] * m.fluxes.r2.value,\n bound = -ParameterBetween(p[3], 0),\n )","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"substitute params into model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"parameter_substitutions = Dict(\n r2bound => 4.0,\n m3bound => 0.1, # lose some mass here\n p[1] => 1.0,\n p[2] => 1.0,\n p[3] => 4.0,\n)\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Tulip.Optimizer,\n)\nm_noparams.fluxes","category":"page"},{"location":"1-parametric-models/#Change-the-parameters-and-re-solve","page":"Parametric constraint-based metabolic models","title":"Change the parameters and re-solve","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"substitute params into model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"parameter_substitutions[m3bound] = 0.0\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Tulip.Optimizer,\n)\nm_noparams.fluxes","category":"page"},{"location":"1-parametric-models/#Quadratic-parameters-also-work","page":"Parametric constraint-based metabolic models","title":"Quadratic parameters also work","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Symbolics.@variables q[1:6]\n\nm.objective = ConstraintTrees.Constraint(\n value = sum(\n rxn.value * rxn.value * qi for (qi, rxn) in zip(collect(q), values(m.fluxes))\n ),\n bound = nothing,\n)\n\nm *= :objective_bound^ConstraintTrees.Constraint(value = m.fluxes.r6.value, bound = 2.0)\n\nparameter_substitutions = merge(parameter_substitutions, Dict(zip(q, fill(1.0, 6))))\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Clarabel.Optimizer,\n sense = Minimal,\n)\nm_noparams.fluxes\n\npqv17 = ParameterLinearValue([1, 2], [2, 1]) * ParameterLinearValue([2], [1]) #src","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"This page was generated using Literate.jl.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"EditURL = \"2-differentiate-enzyme-model.jl\"","category":"page"},{"location":"2-differentiate-enzyme-model/#Differentiating-enzyme-constrained-metabolic-models","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"","category":"section"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"using DifferentiableMetabolism\nusing AbstractFBCModels\nusing Symbolics\nusing ConstraintTrees\nusing COBREXA\nusing Tulip\nusing JSONFBCModels\nimport Downloads: download\nusing CairoMakie\n\n!isfile(\"e_coli_core.json\") &&\n download(\"http://bigg.ucsd.edu/static/models/e_coli_core.json\", \"e_coli_core.json\")\n\ninclude(\"../../test/data_static.jl\")\n\nmodel = load_model(\"e_coli_core.json\")","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"unconstrain glucose","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"rids = [x[\"id\"] for x in model.reactions]\nglc_idx = first(indexin([\"EX_glc__D_e\"], rids))\nmodel.reactions[glc_idx][\"lower_bound\"] = -1000.0","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"constrain PFL to zero","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"pfl_idx = first(indexin([\"PFL\"], rids))\nmodel.reactions[pfl_idx][\"upper_bound\"] = 0.0\n\nkcats = vcat([Symbolics.@variables $x for x in Symbol.(keys(ecoli_core_reaction_kcats))]...)\n\nrid_kcat = Dict(zip(keys(ecoli_core_reaction_kcats), kcats))\n\nparameter_values =\n Dict(kid => ecoli_core_reaction_kcats[rid] * 3.6 for (rid, kid) in rid_kcat) # change unit to k/h\n\nreaction_isozymes = Dict{String,Dict{String,ParameterIsozyme}}() # a mapping from reaction IDs to isozyme IDs to isozyme structs.\nfor rid in AbstractFBCModels.reactions(model)\n grrs = AbstractFBCModels.reaction_gene_association_dnf(model, rid)\n isnothing(grrs) && continue # skip if no grr available\n haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available\n for (i, grr) in enumerate(grrs)\n d = get!(reaction_isozymes, rid, Dict{String,ParameterIsozyme}())\n d[\"isozyme_\"*string(i)] = ParameterIsozyme(\n gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = rid_kcat[rid],\n kcat_reverse = rid_kcat[rid],\n )\n end\nend\n\ngene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses)\n\nSymbolics.@variables capacitylimitation\nparameter_values[capacitylimitation] = 50.0 # mg enzyme/gDW\n\nkm = build_kinetic_model(\n model;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nec_solution, _, _, _ = optimized_constraints_with_parameters(\n km,\n parameter_values;\n objective = km.objective.value,\n optimizer = Tulip.Optimizer,\n modifications = [COBREXA.set_optimizer_attribute(\"IPM_IterationsLimit\", 10_000)],\n)\n\nec_solution","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"This solution contains many inactive reactions","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(abs.(collect(values(ec_solution.fluxes))))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"And also many inactive gene products.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(abs.(collect(values(ec_solution.gene_product_amounts))))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"With theory, you can show that this introduces flux variability into the solution, making it non-unique, and consequently non-differentiable. To fix this, one need to prune the model, to include only the active reactions and genes. This can be differentiated uniquely. Below we build a pruned kinetic model, by removing all the reactions, metabolites, and genes that are not active.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"flux_zero_tol = 1e-6 # these bounds make a real difference!\ngene_zero_tol = 1e-6\n\npruned_reaction_isozymes =\n prune_reaction_isozymes(reaction_isozymes, ec_solution, gene_zero_tol)\n\npruned_gene_product_molar_masses =\n prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol)\n\npkm = build_kinetic_model(\n prune_model(model, ec_solution, flux_zero_tol, gene_zero_tol);\n reaction_isozymes = pruned_reaction_isozymes,\n gene_product_molar_masses = pruned_gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nec_solution2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(\n pkm,\n parameter_values;\n objective = pkm.objective.value,\n optimizer = Tulip.Optimizer,\n modifications = [COBREXA.set_optimizer_attribute(\"IPM_IterationsLimit\", 10_000)],\n)","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Notice, the solution is exactly the same as before, except that all the inactive elements are gone.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"ec_solution2","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"no zero fluxes","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(abs.(collect(values(ec_solution2.fluxes))))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"no zero genes","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(abs.(collect(values(ec_solution2.gene_product_amounts))))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"prune the kcats, leaving only those that are actually used","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"pruned_kcats = [\n begin\n x = first(values(v))\n isnothing(x.kcat_forward) ? x.kcat_reverse : x.kcat_forward\n end for v in values(pruned_reaction_isozymes)\n]\n\nparameters = [capacitylimitation; kcats]\n\nsens, vids = differentiate(\n pkm,\n pkm.objective.value,\n x_vals,\n eq_dual_vals,\n ineq_dual_vals,\n parameter_values,\n parameters;\n scale = true, # unitless sensitivities\n)","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"look at glycolysis and oxidative phosphorylation only","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"subset_ids = [\n r[\"id\"] for r in model.reactions[indexin(string.(keys(pkm.fluxes)), rids)] if\n r[\"subsystem\"] in [\"Glycolysis/Gluconeogenesis\", \"Oxidative Phosphorylation\"]\n]\n\nflux_idxs = findall(x -> string(last(x)) in subset_ids && first(x) == :fluxes, vids)\nflux_ids = last.(vids[flux_idxs])\n\niso_idxs = findall(x -> string(x[2]) in subset_ids && first(x) != :fluxes, vids)\niso_ids = [v[2] for v in vids[iso_idxs]]\n\nparam_idxs = findall(x -> string(x) in subset_ids, parameters)\nparam_ids = parameters[param_idxs]","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Flux sensitivities","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"f, a, hm = heatmap(\n sens[flux_idxs, param_idxs]';\n axis = (\n xlabel = \"kcat\",\n xticks = (1:length(param_ids), string.(param_ids)),\n xticklabelrotation = -pi / 2,\n ylabel = \"Flux\",\n yticks = (1:length(flux_ids), string.(flux_ids)),\n title = \"Flux sensitivities\",\n ),\n)\nColorbar(f[1, 2], hm)\nf","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Isozyme sensitivities. Note, the gene products themselves are not variables in the formulation of the kinetic model. It inherits its structure from COBREXA, where the gene products are derived variables. If you want the sensitivities of the gene products themselves, you just need to multiply the isozyme sensitivity with the subunit stoichiometry of the relevant gene products.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"f, a, hm = heatmap(\n sens[iso_idxs, param_idxs]';\n axis = (\n xlabel = \"kcat\",\n xticks = (1:length(param_ids), string.(param_ids)),\n xticklabelrotation = -pi / 2,\n ylabel = \"Isozyme\",\n yticks = (1:length(iso_ids), string.(iso_ids)),\n title = \"Isozyme sensitivities\",\n ),\n)\nColorbar(f[1, 2], hm)\nf","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"This page was generated using Literate.jl.","category":"page"},{"location":"reference/#Reference","page":"Reference","title":"Reference","text":"","category":"section"},{"location":"reference/#Parameters","page":"Reference","title":"Parameters","text":"","category":"section"},{"location":"reference/#ParameterLinearValue","page":"Reference","title":"ParameterLinearValue","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_linearvalue.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterLinearValue","page":"Reference","title":"DifferentiableMetabolism.ParameterLinearValue","text":"struct ParameterLinearValue <: ConstraintTrees.Value\n\nAn extension of ConstraintTrees.LinearValue where the weights are parameters.\n\nParameterLinearValues can be combined additively and multiplied by real-number constants.\n\nMultiplying two ParameterLinearValues yields a quadratic form (in a ParameterQuadraticValue).\n\nFields\n\nidxs::Vector{Int64}\nweights::Vector{Symbolics.Num}\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterQuadraticValue","page":"Reference","title":"ParameterQuadraticValue","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_quadraticvalue.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterQuadraticValue","page":"Reference","title":"DifferentiableMetabolism.ParameterQuadraticValue","text":"struct ParameterQuadraticValue <: ConstraintTrees.Value\n\nAn extension of ConstraintTrees.QuadraticValue where the weights are parameters.\n\nBehaves similarly to ConstraintTrees.QuadraticValue. Thus, the cleanest way to construct a ParameterQuadraticValue is to multiply two ParameterLinearValues.\n\nFields\n\nidxs::Vector{Tuple{Int64, Int64}}\nweights::Vector{Symbolics.Num}\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterBound","page":"Reference","title":"ParameterBound","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_bound.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterBetween","page":"Reference","title":"DifferentiableMetabolism.ParameterBetween","text":"mutable struct ParameterBetween <: ConstraintTrees.Bound\n\nRepresentation of an \"interval\" bound where the lower and upper bound values are parameters. Since Symbolics.Num is a subtype of Real, the bounds could also be any real number, but they are converted by the constructors to Symbolics.Nums. \n\nFields\n\nlower::Symbolics.Num\nupper::Symbolics.Num\n\n\n\n\n\n","category":"type"},{"location":"reference/#DifferentiableMetabolism.ParameterEqualTo","page":"Reference","title":"DifferentiableMetabolism.ParameterEqualTo","text":"mutable struct ParameterEqualTo <: ConstraintTrees.Bound\n\nRepresentation of an \"equality\" bound, where the bound value is a parameter. Since Symbolics.Num is a subtype of Real, the bound could also be any real number, but it is converted by the constructor to a Symbolics.Num.\n\nFields\n\nequal_to::Symbolics.Num\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterIsozyme","page":"Reference","title":"ParameterIsozyme","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_isozyme.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterIsozyme","page":"Reference","title":"DifferentiableMetabolism.ParameterIsozyme","text":"mutable struct ParameterIsozyme\n\nA parameterized isozyme struct which includes parameters in the kcat_forward, and the kcat_backward. If the reaction does not have a turnover number,nothing can be used. \n\nFields\n\ngene_product_stoichiometry::Dict{String, Float64}\nkcat_forward::Union{Nothing, Symbolics.Num}\nkcat_reverse::Union{Nothing, Symbolics.Num}\n\n\n\n\n\n","category":"type"},{"location":"reference/#Parameterized-models","page":"Reference","title":"Parameterized models","text":"","category":"section"},{"location":"reference/#Kinetic-models","page":"Reference","title":"Kinetic models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/kinetic_model.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.build_kinetic_model-Union{Tuple{AbstractFBCModels.AbstractFBCModel}, Tuple{R}} where R<:Real","page":"Reference","title":"DifferentiableMetabolism.build_kinetic_model","text":"build_kinetic_model(\n model::AbstractFBCModels.AbstractFBCModel;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity\n)\n\n\nEssentially an enzyme constrained metabolic model, but the kcat can be a arbitrary symbolic function.\n\n\n\n\n\n","category":"method"},{"location":"reference/#Solving-models","page":"Reference","title":"Solving models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/solver.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.constraint_matrix_vector-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.constraint_matrix_vector","text":"constraint_matrix_vector(\n eqs,\n m,\n parameters\n) -> Tuple{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseVector{Float64, Int64}}\n\n\nBuilds a matrix representation of bounds.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.optimization_model_with_parameters-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Dict{Symbolics.Num, Float64}}","page":"Reference","title":"DifferentiableMetabolism.optimization_model_with_parameters","text":"optimization_model_with_parameters(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n parameters::Dict{Symbolics.Num, Float64};\n objective,\n optimizer,\n sense\n)\n\n\nConstruct a JuMP model by substituting parameters into the model, m. Set the objective and the optimizer, as well as the sense similar to COBREXA.optimization_model.\n\nConverts all inequality constraints to the form A * x ≤ b.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.optimized_constraints_with_parameters-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Dict{Symbolics.Num, Float64}}","page":"Reference","title":"DifferentiableMetabolism.optimized_constraints_with_parameters","text":"optimized_constraints_with_parameters(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n parameters::Dict{Symbolics.Num, Float64};\n modifications,\n objective,\n optimizer,\n sense\n)\n\n\nSolve a model, m, by forwarding arguments to optimization_model_with_parameters. \n\nOptionally, set optimizer attributes with modifications. If the model does not solve successfully return nothing. Otherwise, return a tuple of the solution tree, and vectors containing the values of the primal variables, the equality constraint dual variables. \n\nThese duals are ordered according to the constraint output of calling equality_constraints and inequality_constraints respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#Pruning-models","page":"Reference","title":"Pruning models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/prune.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.prune_gene_product_molar_masses-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_gene_product_molar_masses","text":"prune_gene_product_molar_masses(\n gene_product_molar_masses,\n ec_solution,\n gene_zero_tol\n) -> Dict\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.prune_model-NTuple{4, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_model","text":"prune_model(\n model,\n ec_solution,\n flux_zero_tol,\n gene_zero_tol\n) -> AbstractFBCModels.CanonicalModel.Model\n\n\nPrune away reactions, metabolites, and genes from a model using ec_solution, which is the result of an enzyme constrained kinetic simulation. Fluxes and gene product concentrations smaller than flux_zero_tol, gene_zero_tol are removed. Metabolites that do not take part in the remaining reactions are also removed.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.prune_reaction_isozymes-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_reaction_isozymes","text":"prune_reaction_isozymes(\n reaction_isozymes,\n ec_solution,\n gene_zero_tol\n) -> Dict\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#Differentiation","page":"Reference","title":"Differentiation","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/differentiate.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.differentiate-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, ConstraintTrees.Value, Vector{Float64}, Vector{Float64}, Vector{Float64}, Dict{Symbolics.Num, Float64}, Vector{Symbolics.Num}}","page":"Reference","title":"DifferentiableMetabolism.differentiate","text":"differentiate(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n objective::ConstraintTrees.Value,\n x_vals::Vector{Float64},\n eq_dual_vals::Vector{Float64},\n ineq_dual_vals::Vector{Float64},\n parameter_values::Dict{Symbolics.Num, Float64},\n parameters::Vector{Symbolics.Num};\n scale\n) -> Tuple{Any, Any}\n\n\nDifferentiate a model m with respect to parameters which take on values parameter_values in the optimal solution with respect to the objective. The primal variables x_vals, and the dual variable values eq_dual_vals and ineq_dual_vals need to be supplied. \n\nInternally, primal variables with value abs(x) ≤ primal_zero_tol are removed from the computation, and their sensitivities are not calculated. \n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.findall_indeps_qr-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.findall_indeps_qr","text":"findall_indeps_qr(A) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.variable_order-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.variable_order","text":"variable_order(m) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#Internals","page":"Reference","title":"Internals","text":"","category":"section"},{"location":"reference/#Constraint-tree-extensions","page":"Reference","title":"Constraint tree extensions","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/constraint_trees.jl\"]","category":"page"},{"location":"reference/#Optimization-problem-builders","page":"Reference","title":"Optimization problem builders","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/get_constraints.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.constrained-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.constrained","text":"constrained(x) -> Bool\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.equality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}}","page":"Reference","title":"DifferentiableMetabolism.equality_constraints","text":"equality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint}\n) -> Vector{Tuple{Union{ConstraintTrees.LinearValue, ParameterLinearValue}, Symbolics.Num}}\n\n\nReturn all the equality constraints of m as a tuple ({Parameter}LinearValue, value) representing {P}LV == value for each entry.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_equality_constraints-Tuple{ConstraintTrees.Constraint, Any}","page":"Reference","title":"DifferentiableMetabolism.get_equality_constraints","text":"get_equality_constraints(\n c::ConstraintTrees.Constraint,\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_equality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Any}","page":"Reference","title":"DifferentiableMetabolism.get_equality_constraints","text":"get_equality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_inequality_constraints-Tuple{ConstraintTrees.Constraint, Any}","page":"Reference","title":"DifferentiableMetabolism.get_inequality_constraints","text":"get_inequality_constraints(\n c::ConstraintTrees.Constraint,\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_inequality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Any}","page":"Reference","title":"DifferentiableMetabolism.get_inequality_constraints","text":"get_inequality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.inequality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}}","page":"Reference","title":"DifferentiableMetabolism.inequality_constraints","text":"inequality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint}\n) -> Vector\n\n\nReturn all the inequality constraints of m as a tuple of bounds converted to the form ({Parameter}LinearValue, upper) representing {P}LV ≤ upper for each entry. \n\n\n\n\n\n","category":"method"},{"location":"reference/#Symbolics-extensions","page":"Reference","title":"Symbolics extensions","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/symbolics.jl\"]","category":"page"},{"location":"#DifferentiableMetabolism.jl","page":"README","title":"DifferentiableMetabolism.jl","text":"","category":"section"},{"location":"","page":"README","title":"README","text":"Documentation for DifferentiableMetabolism.","category":"page"},{"location":"","page":"README","title":"README","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/DifferentiableMetabolism.jl\"]","category":"page"}] +[{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"EditURL = \"3-parameter-estimation.jl\"","category":"page"},{"location":"3-parameter-estimation/#Parameter-estimation-using-proteomics-and-flux-data","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"","category":"section"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"using DifferentiableMetabolism\nusing AbstractFBCModels\nusing FastDifferentiation\nusing ConstraintTrees\nusing COBREXA\nusing Clarabel\nusing Tulip\nusing JSONFBCModels\nusing CairoMakie\n\ninclude(\"../../test/simple_model.jl\") #hide","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"prune model for brevity","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"delete!(model.reactions, \"r5\")\ndelete!(model.genes, \"g4\")\ndelete!(model.genes, \"g5\")\ndelete!(model.genes, \"g3\")\nmodel.reactions[\"r4\"].gene_association_dnf = [[\"g2\"]]\nmodel.reactions[\"r1\"].lower_bound = -1000.0\nmodel.reactions[\"r2\"].lower_bound = -1000.0","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"now models looks like this","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"(Image: simple_model)","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"@variables r3 r4\n\nreaction_isozymes = Dict(\n \"r3\" => Dict(\n \"isozyme1\" => ParameterIsozyme(\n gene_product_stoichiometry = Dict(\"g1\" => 1.0), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = r3,\n kcat_reverse = nothing,\n ),\n ),\n \"r4\" => Dict(\n \"isozyme1\" => ParameterIsozyme(\n gene_product_stoichiometry = Dict(\"g2\" => 1.0), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = r4,\n kcat_reverse = nothing,\n ),\n ),\n)\n\ngene_product_molar_masses = Dict(\"g1\" => 20.0, \"g2\" => 10.0)\n\n@variables capacitylimitation\n\ntrue_parameter_values = Dict(:capacitylimitation => 50.0, :r3 => 2.0, :r4 => 3.0)\n\nkm = build_kinetic_model(\n model;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nsol, _, _, _ = optimized_constraints_with_parameters(\n km,\n true_parameter_values;\n objective = km.objective.value,\n optimizer = Tulip.Optimizer,\n)\n\nsol.fluxes","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"create a loss function","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"measured = [\n sol.fluxes.r1,\n sol.fluxes.r3,\n sol.isozyme_forward_amounts.r3.isozyme1,\n sol.isozyme_forward_amounts.r4.isozyme1,\n]\n\nkm *=\n :loss^ConstraintTrees.Constraint(;\n value = 0.5 * (\n ConstraintTrees.squared(km.fluxes.r1.value - measured[1]) +\n ConstraintTrees.squared(km.fluxes.r3.value - measured[2]) +\n ConstraintTrees.squared(\n km.isozyme_forward_amounts.r3.isozyme1.value - measured[3],\n ) +\n ConstraintTrees.squared(\n km.isozyme_forward_amounts.r4.isozyme1.value - measured[4],\n )\n ),\n bound = nothing,\n )\n\nestimated_parameters = Dict(:capacitylimitation => 50.0, :r3 => 5.0, :r4 => 1.0) # initial values\nη = 0.1 # learning rate\n\nlosses = Float64[]\n\nkmKKT, vids = differentiate_prepare_kkt(\n km,\n km.loss.value,\n [:r3, :r4, :capacitylimitation],\n)\n\nfor k = 1:150\n\n sol2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(\n km,\n estimated_parameters;\n objective = km.loss.value,\n optimizer = Clarabel.Optimizer,\n sense = COBREXA.Minimal,\n modifications = [COBREXA.silence],\n )\n\n push!(losses, sol2.loss)\n\n sens = differentiate_solution(\n kmKKT,\n x_vals,\n eq_dual_vals,\n ineq_dual_vals,\n estimated_parameters,\n )\n measured_idxs = [1, 3, 12, 11]\n\n x = [\n sol2.fluxes.r1,\n sol2.fluxes.r3,\n sol2.isozyme_forward_amounts.r3.isozyme1,\n sol2.isozyme_forward_amounts.r4.isozyme1,\n ]\n\n dL_dx = x - measured # derivative of loss function with respect to optimization variables\n dL_dkcats = sens[measured_idxs, :]' * dL_dx # derivative of loss function with respect to parameters\n\n estimated_parameters[:r3] -= η * dL_dkcats[1]\n estimated_parameters[:r4] -= η * dL_dkcats[2]\nend\n\nlines(losses; axis = (xlabel = \"Iterations\", ylabel = \"L2 loss\"))","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"estimated_parameters","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"true_parameter_values","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"","category":"page"},{"location":"3-parameter-estimation/","page":"Parameter estimation using proteomics and flux data","title":"Parameter estimation using proteomics and flux data","text":"This page was generated using Literate.jl.","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"EditURL = \"1-parametric-models.jl\"","category":"page"},{"location":"1-parametric-models/#Parametric-constraint-based-metabolic-models","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"using DifferentiableMetabolism\n\nusing FastDifferentiation\nusing ConstraintTrees\nusing COBREXA\nusing Tulip\nusing Clarabel","category":"page"},{"location":"1-parametric-models/#Load-and-solve-a-simple-model","page":"Parametric constraint-based metabolic models","title":"Load and solve a simple model","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"The code used to construct the model is located in test/simple_model.jl, but it is not shown here for brevity. Below is a visualization of the model.","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"include(\"../../test/simple_model.jl\"); #hide\nnothing #hide","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"(Image: simple_model)","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Build a basic ConstraintTree model without parameters","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"m = COBREXA.flux_balance_constraints(model)","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Solve normally","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"base_model =\n COBREXA.optimized_values(m; optimizer = Tulip.Optimizer, objective = m.objective.value)\nbase_model.fluxes","category":"page"},{"location":"1-parametric-models/#Add-parameters-to-the-model","page":"Parametric constraint-based metabolic models","title":"Add parameters to the model","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"Make bound of r2 and mass balance of m3 parameters","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"@variables r2bound m3bound\n\nm.fluxes.r2 =\n ConstraintTrees.Constraint(m.fluxes.r2.value, -2 * ParameterBetween(r2bound, 0))\n\nm.flux_stoichiometry.m3 =\n ConstraintTrees.Constraint(m.flux_stoichiometry.m3.value, ParameterEqualTo(m3bound) / 2)","category":"page"},{"location":"1-parametric-models/#add-parametric-constraints","page":"Parametric constraint-based metabolic models","title":"add parametric constraints","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"p = make_variables(:p, 4)\n\nm *=\n :linparam^ConstraintTrees.Constraint(\n value = p[1] * m.fluxes.r1.value + p[2] * m.fluxes.r2.value,\n bound = -ParameterBetween(p[3], 0),\n )","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"substitute params into model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"parameter_substitutions = Dict(\n :r2bound => 4.0,\n :m3bound => 0.1, # lose some mass here\n :p1 => 1.0,\n :p2 => 1.0,\n :p3 => 4.0,\n)\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Tulip.Optimizer,\n)\nm_noparams.fluxes","category":"page"},{"location":"1-parametric-models/#Change-the-parameters-and-re-solve","page":"Parametric constraint-based metabolic models","title":"Change the parameters and re-solve","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"substitute params into model","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"parameter_substitutions[:m3bound] = 0.0\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Tulip.Optimizer,\n)\nm_noparams.fluxes","category":"page"},{"location":"1-parametric-models/#Quadratic-parameters-also-work","page":"Parametric constraint-based metabolic models","title":"Quadratic parameters also work","text":"","category":"section"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"q = make_variables(:q, 6)\n\nm.objective = ConstraintTrees.Constraint(\n value = sum(\n rxn.value * rxn.value * qi for (qi, rxn) in zip(collect(q), values(m.fluxes))\n ),\n bound = nothing,\n)\n\nm *= :objective_bound^ConstraintTrees.Constraint(value = m.fluxes.r6.value, bound = 2.0)\n\nparameter_substitutions =\n merge(parameter_substitutions, Dict(v.node_value => 1.0 for v in q))\n\nm_noparams, _, _, _ = optimized_constraints_with_parameters(\n m,\n parameter_substitutions;\n objective = m.objective.value,\n optimizer = Clarabel.Optimizer,\n sense = Minimal,\n)\nm_noparams.fluxes","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"","category":"page"},{"location":"1-parametric-models/","page":"Parametric constraint-based metabolic models","title":"Parametric constraint-based metabolic models","text":"This page was generated using Literate.jl.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"EditURL = \"2-differentiate-enzyme-model.jl\"","category":"page"},{"location":"2-differentiate-enzyme-model/#Differentiating-enzyme-constrained-metabolic-models","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"","category":"section"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"using DifferentiableMetabolism\nusing AbstractFBCModels\nusing ConstraintTrees\nusing COBREXA\nusing Tulip\nusing JSONFBCModels\nimport Downloads: download\nusing CairoMakie\nusing FastDifferentiation\n\n!isfile(\"e_coli_core.json\") &&\n download(\"http://bigg.ucsd.edu/static/models/e_coli_core.json\", \"e_coli_core.json\")\n\ninclude(\"../../test/data_static.jl\")\n\nmodel = load_model(\"e_coli_core.json\")","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"unconstrain glucose","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"rids = [x[\"id\"] for x in model.reactions]\nglc_idx = first(indexin([\"EX_glc__D_e\"], rids))\nmodel.reactions[glc_idx][\"lower_bound\"] = -1000.0","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"constrain PFL to zero","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"pfl_idx = first(indexin([\"PFL\"], rids))\nmodel.reactions[pfl_idx][\"upper_bound\"] = 0.0\n\nrid_kcat = Dict(k => FastDifferentiation.Node(Symbol(k)) for (k,_) in ecoli_core_reaction_kcats)\nkcats = Symbol.(keys(ecoli_core_reaction_kcats))\n\nparameter_values = Dict{Symbol, Float64}()\n\nreaction_isozymes = Dict{String,Dict{String,ParameterIsozyme}}() # a mapping from reaction IDs to isozyme IDs to isozyme structs.\nfor rid in AbstractFBCModels.reactions(model)\n grrs = AbstractFBCModels.reaction_gene_association_dnf(model, rid)\n isnothing(grrs) && continue # skip if no grr available\n haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available\n for (i, grr) in enumerate(grrs)\n\n kcat = ecoli_core_reaction_kcats[rid] * 3.6 # change unit to k/h\n parameter_values[Symbol(rid)] = kcat\n\n d = get!(reaction_isozymes, rid, Dict{String,ParameterIsozyme}())\n d[\"isozyme_$i\"] = ParameterIsozyme(\n gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes\n kcat_forward = rid_kcat[rid],\n kcat_reverse = rid_kcat[rid],\n )\n end\nend\n\ngene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses)\n\n@variables capacitylimitation\nparameter_values[:capacitylimitation] = 50.0 # mg enzyme/gDW\n\nkm = build_kinetic_model(\n model;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nec_solution, _, _, _ = optimized_constraints_with_parameters(\n km,\n parameter_values;\n objective = km.objective.value,\n optimizer = Tulip.Optimizer,\n modifications = [COBREXA.set_optimizer_attribute(\"IPM_IterationsLimit\", 10_000)],\n)\n\nec_solution","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"This solution contains many inactive reactions","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"And also many inactive gene products.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(collect(ec_solution.gene_product_amounts), by=last)","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"With theory, you can show that this introduces flux variability into the solution, making it non-unique, and consequently non-differentiable. To fix this, one need to prune the model, to include only the active reactions and genes. This can be differentiated uniquely. Below we build a pruned kinetic model, by removing all the reactions, metabolites, and genes that are not active.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"flux_zero_tol = 1e-6 # these bounds make a real difference!\ngene_zero_tol = 1e-6\n\npruned_reaction_isozymes =\n prune_reaction_isozymes(reaction_isozymes, ec_solution, flux_zero_tol)\n\npruned_gene_product_molar_masses =\n prune_gene_product_molar_masses(gene_product_molar_masses, ec_solution, gene_zero_tol)\n\npkm = build_kinetic_model(\n prune_model(model, ec_solution, flux_zero_tol, gene_zero_tol);\n reaction_isozymes = pruned_reaction_isozymes,\n gene_product_molar_masses = pruned_gene_product_molar_masses,\n capacity = capacitylimitation,\n)\n\nec_solution2, x_vals, eq_dual_vals, ineq_dual_vals = optimized_constraints_with_parameters(\n pkm,\n parameter_values;\n objective = pkm.objective.value,\n optimizer = Tulip.Optimizer,\n modifications = [COBREXA.set_optimizer_attribute(\"IPM_IterationsLimit\", 10_000)],\n)","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Notice, the solution is exactly the same as before, except that all the inactive elements are gone.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"ec_solution2","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"no zero fluxes","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(collect(ec_solution.fluxes), by=ComposedFunction(abs, last))","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"no zero genes","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"sort(abs.(collect(values(ec_solution2.gene_product_amounts))))\n\n\n\n\nparameters = [:capacitylimitation; kcats]\n\npkm_kkt, vids = differentiate_prepare_kkt(\n pkm,\n pkm.objective.value,\n parameters\n)\n\nsens = differentiate_solution(\n pkm_kkt,\n x_vals,\n eq_dual_vals,\n ineq_dual_vals,\n parameter_values,\n scale = true, # unitless sensitivities\n)","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"look at glycolysis and oxidative phosphorylation only","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"subset_ids = [\n r[\"id\"] for r in model.reactions[indexin(string.(keys(pkm.fluxes)), rids)] if\n r[\"subsystem\"] in [\"Glycolysis/Gluconeogenesis\", \"Oxidative Phosphorylation\"]\n]\n\nflux_idxs = findall(x -> string(last(x)) in subset_ids && first(x) == :fluxes, vids)\nflux_ids = last.(vids[flux_idxs])\n\niso_idxs = findall(x -> string(x[2]) in subset_ids && first(x) != :fluxes, vids)\niso_ids = [v[2] for v in vids[iso_idxs]]\n\nparam_idxs = findall(x -> string(x) in subset_ids, parameters)\nparam_ids = parameters[param_idxs]","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Flux sensitivities","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"f, a, hm = heatmap(\n sens[flux_idxs, param_idxs]';\n axis = (\n xlabel = \"kcat\",\n xticks = (1:length(param_ids), string.(param_ids)),\n xticklabelrotation = -pi / 2,\n ylabel = \"Flux\",\n yticks = (1:length(flux_ids), string.(flux_ids)),\n title = \"Flux sensitivities\",\n ),\n)\nColorbar(f[1, 2], hm)\nf","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"Isozyme sensitivities. Note, the gene products themselves are not variables in the formulation of the kinetic model. It inherits its structure from COBREXA, where the gene products are derived variables. If you want the sensitivities of the gene products themselves, you just need to multiply the isozyme sensitivity with the subunit stoichiometry of the relevant gene products.","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"f, a, hm = heatmap(\n sens[iso_idxs, param_idxs]';\n axis = (\n xlabel = \"kcat\",\n xticks = (1:length(param_ids), string.(param_ids)),\n xticklabelrotation = -pi / 2,\n ylabel = \"Isozyme\",\n yticks = (1:length(iso_ids), string.(iso_ids)),\n title = \"Isozyme sensitivities\",\n ),\n)\nColorbar(f[1, 2], hm)\nf","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"","category":"page"},{"location":"2-differentiate-enzyme-model/","page":"Differentiating enzyme constrained metabolic models","title":"Differentiating enzyme constrained metabolic models","text":"This page was generated using Literate.jl.","category":"page"},{"location":"reference/#Reference","page":"Reference","title":"Reference","text":"","category":"section"},{"location":"reference/#Parameters","page":"Reference","title":"Parameters","text":"","category":"section"},{"location":"reference/#ParameterLinearValue","page":"Reference","title":"ParameterLinearValue","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_linearvalue.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterLinearValue","page":"Reference","title":"DifferentiableMetabolism.ParameterLinearValue","text":"struct ParameterLinearValue <: ConstraintTrees.Value\n\nAn extension of ConstraintTrees.LinearValue where the weights are parameters.\n\nParameterLinearValues can be combined additively and multiplied by real-number constants.\n\nMultiplying two ParameterLinearValues yields a quadratic form (in a ParameterQuadraticValue).\n\nFields\n\nidxs::Vector{Int64}\nweights::Vector{FastDifferentiation.Node}\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterQuadraticValue","page":"Reference","title":"ParameterQuadraticValue","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_quadraticvalue.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterQuadraticValue","page":"Reference","title":"DifferentiableMetabolism.ParameterQuadraticValue","text":"struct ParameterQuadraticValue <: ConstraintTrees.Value\n\nAn extension of ConstraintTrees.QuadraticValue where the weights are parameters.\n\nBehaves similarly to ConstraintTrees.QuadraticValue. Thus, the cleanest way to construct a ParameterQuadraticValue is to multiply two ParameterLinearValues.\n\nFields\n\nidxs::Vector{Tuple{Int64, Int64}}\nweights::Vector{FastDifferentiation.Node}\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterBound","page":"Reference","title":"ParameterBound","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_bound.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterBetween","page":"Reference","title":"DifferentiableMetabolism.ParameterBetween","text":"mutable struct ParameterBetween <: ConstraintTrees.Bound\n\nRepresentation of an \"interval\" bound where the lower and upper bound values are parameters. Since Expression is a subtype of Real, the bounds could also be any real number, but they are converted by the constructors to Expressions. \n\nFields\n\nlower::FastDifferentiation.Node\nupper::FastDifferentiation.Node\n\n\n\n\n\n","category":"type"},{"location":"reference/#DifferentiableMetabolism.ParameterEqualTo","page":"Reference","title":"DifferentiableMetabolism.ParameterEqualTo","text":"mutable struct ParameterEqualTo <: ConstraintTrees.Bound\n\nRepresentation of an \"equality\" bound, where the bound value is a parameter. Since Expression is a subtype of Real, the bound could also be any real number, but it is converted by the constructor to a Expression.\n\nFields\n\nequal_to::FastDifferentiation.Node\n\n\n\n\n\n","category":"type"},{"location":"reference/#ParameterIsozyme","page":"Reference","title":"ParameterIsozyme","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/parameter_isozyme.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.ParameterIsozyme","page":"Reference","title":"DifferentiableMetabolism.ParameterIsozyme","text":"mutable struct ParameterIsozyme\n\nA parameterized isozyme struct which includes parameters in the kcat_forward, and the kcat_backward. If the reaction does not have a turnover number,nothing can be used. \n\nFields\n\ngene_product_stoichiometry::Dict{String, Float64}\nkcat_forward::Union{Nothing, FastDifferentiation.Node}\nkcat_reverse::Union{Nothing, FastDifferentiation.Node}\n\n\n\n\n\n","category":"type"},{"location":"reference/#Parameterized-models","page":"Reference","title":"Parameterized models","text":"","category":"section"},{"location":"reference/#Kinetic-models","page":"Reference","title":"Kinetic models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/kinetic_model.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.build_kinetic_model-Union{Tuple{AbstractFBCModels.AbstractFBCModel}, Tuple{R}} where R<:Real","page":"Reference","title":"DifferentiableMetabolism.build_kinetic_model","text":"build_kinetic_model(\n model::AbstractFBCModels.AbstractFBCModel;\n reaction_isozymes,\n gene_product_molar_masses,\n capacity\n)\n\n\nEssentially an enzyme constrained metabolic model, but the kcat can be a arbitrary symbolic function.\n\n\n\n\n\n","category":"method"},{"location":"reference/#Solving-models","page":"Reference","title":"Solving models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/solver.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.constraint_matrix_vector-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.constraint_matrix_vector","text":"constraint_matrix_vector(\n eqs,\n m,\n parameters\n) -> Tuple{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseVector{Float64, Int64}}\n\n\nBuilds a matrix representation of bounds.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.optimization_model_with_parameters-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Dict{Symbol, Float64}}","page":"Reference","title":"DifferentiableMetabolism.optimization_model_with_parameters","text":"optimization_model_with_parameters(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n parameters::Dict{Symbol, Float64};\n objective,\n optimizer,\n sense\n)\n\n\nConstruct a JuMP model by substituting parameters into the model, m. Set the objective and the optimizer, as well as the sense similar to COBREXA.optimization_model.\n\nConverts all inequality constraints to the form A * x ≤ b.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.optimized_constraints_with_parameters-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Dict{Symbol, Float64}}","page":"Reference","title":"DifferentiableMetabolism.optimized_constraints_with_parameters","text":"optimized_constraints_with_parameters(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n parameters::Dict{Symbol, Float64};\n modifications,\n objective,\n optimizer,\n sense\n)\n\n\nSolve a model, m, by forwarding arguments to optimization_model_with_parameters. \n\nOptionally, set optimizer attributes with modifications. If the model does not solve successfully return nothing. Otherwise, return a tuple of the solution tree, and vectors containing the values of the primal variables, the equality constraint dual variables. \n\nThese duals are ordered according to the constraint output of calling equality_constraints and inequality_constraints respectively.\n\n\n\n\n\n","category":"method"},{"location":"reference/#Pruning-models","page":"Reference","title":"Pruning models","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/prune.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.prune_gene_product_molar_masses-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_gene_product_molar_masses","text":"prune_gene_product_molar_masses(\n gene_product_molar_masses,\n ec_solution,\n gene_zero_tol\n) -> Dict\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.prune_model-NTuple{4, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_model","text":"prune_model(\n model,\n ec_solution,\n flux_zero_tol,\n gene_zero_tol\n) -> AbstractFBCModels.CanonicalModel.Model\n\n\nPrune away reactions, metabolites, and genes from a model using ec_solution, which is the result of an enzyme constrained kinetic simulation. Fluxes and gene product concentrations smaller than flux_zero_tol, gene_zero_tol are removed. Metabolites that do not take part in the remaining reactions are also removed.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.prune_reaction_isozymes-Tuple{Any, Any, Any}","page":"Reference","title":"DifferentiableMetabolism.prune_reaction_isozymes","text":"prune_reaction_isozymes(\n reaction_isozymes,\n ec_solution,\n gene_zero_tol\n) -> Dict\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#Differentiation","page":"Reference","title":"Differentiation","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/differentiate.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.differentiate_prepare_kkt-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, ConstraintTrees.Value, Vector{Symbol}}","page":"Reference","title":"DifferentiableMetabolism.differentiate_prepare_kkt","text":"differentiate_prepare_kkt(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n objective::ConstraintTrees.Value,\n parameters::Vector{Symbol}\n) -> Tuple{Tuple{SparseArrays.SparseMatrixCSC{FastDifferentiation.Node, Int64}, SparseArrays.SparseMatrixCSC{FastDifferentiation.Node, Int64}, Any, Vector{FastDifferentiation.Node}, Vector{FastDifferentiation.Node}, Vector{Symbol}}, Any}\n\n\nPrepare a model m with objective for differentiation with respect to parameters.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.findall_indeps_qr-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.findall_indeps_qr","text":"findall_indeps_qr(A) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.variable_order-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.variable_order","text":"variable_order(m) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#Internals","page":"Reference","title":"Internals","text":"","category":"section"},{"location":"reference/#Constraint-tree-extensions","page":"Reference","title":"Constraint tree extensions","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/constraint_trees.jl\"]","category":"page"},{"location":"reference/#Optimization-problem-builders","page":"Reference","title":"Optimization problem builders","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/get_constraints.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.constrained-Tuple{Any}","page":"Reference","title":"DifferentiableMetabolism.constrained","text":"constrained(x) -> Bool\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.equality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}}","page":"Reference","title":"DifferentiableMetabolism.equality_constraints","text":"equality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint}\n) -> Vector{Tuple{Union{ConstraintTrees.LinearValue, ParameterLinearValue}, FastDifferentiation.Node}}\n\n\nReturn all the equality constraints of m as a tuple ({Parameter}LinearValue, value) representing {P}LV == value for each entry.\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_equality_constraints-Tuple{ConstraintTrees.Constraint, Any}","page":"Reference","title":"DifferentiableMetabolism.get_equality_constraints","text":"get_equality_constraints(\n c::ConstraintTrees.Constraint,\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_equality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Any}","page":"Reference","title":"DifferentiableMetabolism.get_equality_constraints","text":"get_equality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_inequality_constraints-Tuple{ConstraintTrees.Constraint, Any}","page":"Reference","title":"DifferentiableMetabolism.get_inequality_constraints","text":"get_inequality_constraints(\n c::ConstraintTrees.Constraint,\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.get_inequality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}, Any}","page":"Reference","title":"DifferentiableMetabolism.get_inequality_constraints","text":"get_inequality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint},\n sink\n) -> Any\n\n\nTODO\n\n\n\n\n\n","category":"method"},{"location":"reference/#DifferentiableMetabolism.inequality_constraints-Tuple{ConstraintTrees.Tree{ConstraintTrees.Constraint}}","page":"Reference","title":"DifferentiableMetabolism.inequality_constraints","text":"inequality_constraints(\n m::ConstraintTrees.Tree{ConstraintTrees.Constraint}\n) -> Vector\n\n\nReturn all the inequality constraints of m as a tuple of bounds converted to the form ({Parameter}LinearValue, upper) representing {P}LV ≤ upper for each entry. \n\n\n\n\n\n","category":"method"},{"location":"reference/#Expression-evaluation-utilities","page":"Reference","title":"Expression evaluation utilities","text":"","category":"section"},{"location":"reference/","page":"Reference","title":"Reference","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/substitute.jl\"]","category":"page"},{"location":"reference/#DifferentiableMetabolism.substitute-Tuple{FastDifferentiation.Node, Any}","page":"Reference","title":"DifferentiableMetabolism.substitute","text":"substitute(x::FastDifferentiation.Node, lookup) -> Any\n\n\nStraightforward recursive evaluator for Expressions.\n\n\n\n\n\n","category":"method"},{"location":"#DifferentiableMetabolism.jl","page":"README","title":"DifferentiableMetabolism.jl","text":"","category":"section"},{"location":"","page":"README","title":"README","text":"Documentation for DifferentiableMetabolism.","category":"page"},{"location":"","page":"README","title":"README","text":"Modules = [DifferentiableMetabolism]\nPages = [\"src/DifferentiableMetabolism.jl\"]","category":"page"}] }