From 42e378bc0bb893df6e7676406dd89b5790cbde00 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Sat, 12 Oct 2024 23:55:22 -0700 Subject: [PATCH 1/2] Add parameters to HIRES --- toolbox/+otp/+hires/+presets/Canonical.m | 3 +- toolbox/+otp/+hires/HIRESParameters.m | 32 ++++++++++++ toolbox/+otp/+hires/HIRESProblem.m | 22 ++++++-- toolbox/+otp/+hires/f.m | 51 ++++++++++++++----- toolbox/+otp/+hires/jacobian.m | 28 ++++------ .../+hires/jacobianAdjointVectorProduct.m | 34 +++++-------- toolbox/+otp/+hires/jacobianVectorProduct.m | 45 ++++++++-------- 7 files changed, 136 insertions(+), 79 deletions(-) diff --git a/toolbox/+otp/+hires/+presets/Canonical.m b/toolbox/+otp/+hires/+presets/Canonical.m index 5fc99b07..c4d1b49a 100644 --- a/toolbox/+otp/+hires/+presets/Canonical.m +++ b/toolbox/+otp/+hires/+presets/Canonical.m @@ -3,7 +3,8 @@ function obj = Canonical(varargin) tspan = [0, 321.8122]; y0 = [1; 0; 0; 0; 0; 0; 0; 0.0057]; - params = otp.hires.HIRESParameters(varargin{:}); + params = otp.hires.HIRESParameters('K1', 1.71, 'K2', 0.43, 'K3', 8.32, 'K4', 0.69, 'K5', 0.035, ... + 'K6', 8.32, 'KPlus', 280, 'KMinus', 0.69, 'KStar', 0.69, 'OKS', 0.0007, varargin{:}); obj = obj@otp.hires.HIRESProblem(tspan, y0, params); end end diff --git a/toolbox/+otp/+hires/HIRESParameters.m b/toolbox/+otp/+hires/HIRESParameters.m index e4fd66bc..5a8de173 100644 --- a/toolbox/+otp/+hires/HIRESParameters.m +++ b/toolbox/+otp/+hires/HIRESParameters.m @@ -1,6 +1,38 @@ classdef HIRESParameters < otp.Parameters % Parameters for the HIRES problem. + properties + % The stiffness parameter $ε$ for the "cusp catastrophe" model. + K1 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + K2 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + K3 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + K4 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + K5 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + K6 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + KPlus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + KMinus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + KStar %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + + % The diffusion coefficient $σ$ for all three variables. + OKS %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} + end + methods function obj = HIRESParameters(varargin) % Create a HIRES parameters object. diff --git a/toolbox/+otp/+hires/HIRESProblem.m b/toolbox/+otp/+hires/HIRESProblem.m index 1912a8dc..194592db 100644 --- a/toolbox/+otp/+hires/HIRESProblem.m +++ b/toolbox/+otp/+hires/HIRESProblem.m @@ -9,10 +9,24 @@ methods (Access = protected) function onSettingsChanged(obj) - obj.RHS = otp.RHS(@otp.hires.f, ... - 'Jacobian', @otp.hires.jacobian, ... - 'JacobianVectorProduct', @otp.hires.jacobianVectorProduct, ... - 'JacobianAdjointVectorProduct', @otp.hires.jacobianAdjointVectorProduct); + k1 = obj.Parameters.K1; + k2 = obj.Parameters.K2; + k3 = obj.Parameters.K3; + k4 = obj.Parameters.K4; + k5 = obj.Parameters.K5; + k6 = obj.Parameters.K6; + kPlus = obj.Parameters.KPlus; + kMinus = obj.Parameters.KMinus; + kStar = obj.Parameters.KStar; + oks = obj.Parameters.OKS; + + obj.RHS = otp.RHS(@(t, y) otp.hires.f(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... + 'Jacobian', @(t, y) otp.hires.jacobian(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... + 'JacobianVectorProduct', @(t, y, v) ... + otp.hires.jacobianVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... + 'JacobianAdjointVectorProduct', @(t, y, v) ... + otp.hires.jacobianAdjointVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... + 'Vectorized', 'on'); end function fig = internalPlot(obj, t, y, varargin) diff --git a/toolbox/+otp/+hires/f.m b/toolbox/+otp/+hires/f.m index 6caf7f2f..c7a275e0 100644 --- a/toolbox/+otp/+hires/f.m +++ b/toolbox/+otp/+hires/f.m @@ -1,13 +1,38 @@ -function yPrime = f(~, y) - -yPrime = [ - -1.71 * y(1) + 0.43 * y(2) + 8.32 * y(3) + 0.0007; - 1.71 * y(1) - 8.75 * y(2); - -10.03 * y(3) + 0.43 * y(4) + 0.035 * y(5); - 8.32 * y(2) + 1.71 * y(3) - 1.12 * y(4); - -1.745 * y(5) + 0.43 * y(6) + 0.43 * y(7); - -280 * y(6) * y(8) + 0.69 * y(4) + 1.71 * y(5) - 0.43 * y(6) + 0.69 * y(7); - 280 * y(6) * y(8) - 1.81 * y(7); - -280 * y(6) * y(8) + 1.81 * y(7)]; - -end \ No newline at end of file +function yPrime = f(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks) + +k1Pr = k1 * y(1, :); +k2Pfr = k2 * y(2, :); +k3Pfr = k3 * y(2, :); +k1PrX = k1 * y(3, :); +k6PrX = k6 * y(3, :); +k2PfrX = k2 * y(4, :); +k4PfrX = k4 * y(4, :); +k1PrXPrime = k1 * y(5, :); +k5PrXPrime = k5 * y(5, :); +k2PfrXPrime = k2 * y(6, :); +k2PfrXPrimeE = k2 * y(7, :); +kMinusPfrXPrimeE = kMinus * y(7, :); +kPlusPfrXPrimeE = kPlus * y(6, :) * y(8, :); +yPrime7 = -k2PfrXPrimeE - (kMinus + kStar) * y(7, :) + kPlusPfrXPrimeE; + +yPrime = [ ... + -k1Pr + k2Pfr + k6PrX + oks; + k1Pr - k2Pfr - k3Pfr; + -k1PrX - k6PrX + k2PfrX + k5PrXPrime; + k3Pfr + k1PrX - k2PfrX - k4PfrX; + -k1PrXPrime - k5PrXPrime + k2PfrXPrime + k2PfrXPrimeE; + k4PfrX + k1PrXPrime - k2PfrXPrime + kMinusPfrXPrimeE - kPlusPfrXPrimeE; + yPrime7; + -yPrime7]; + +% yPrime = [ +% -1.71 * y(1) + 0.43 * y(2) + 8.32 * y(3) + 0.0007; +% 1.71 * y(1) - 8.75 * y(2); +% -10.03 * y(3) + 0.43 * y(4) + 0.035 * y(5); +% 8.32 * y(2) + 1.71 * y(3) - 1.12 * y(4); +% -1.745 * y(5) + 0.43 * y(6) + 0.43 * y(7); +% -280 * y(6) * y(8) + 0.69 * y(4) + 1.71 * y(5) - 0.43 * y(6) + 0.69 * y(7); +% 280 * y(6) * y(8) - 1.81 * y(7); +% -280 * y(6) * y(8) + 1.81 * y(7)]; + +end diff --git a/toolbox/+otp/+hires/jacobian.m b/toolbox/+otp/+hires/jacobian.m index 429c1e02..b9c69aa3 100644 --- a/toolbox/+otp/+hires/jacobian.m +++ b/toolbox/+otp/+hires/jacobian.m @@ -1,23 +1,13 @@ -function J = jacobian(~, y) +function jac = jacobian(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~) -y6 = y(6); -y8 = y(8); +kPlusPfrXPrime = kPlus * y(6); +kPlusE = kPlus * y(8); +kSum = k2 + kMinus + kStar; -J = sparse( ... - [1,2,1,2,4,1,3,4,3,4,6,3,5,6,5,6,7,8,5,6,7,8,6,7,8], ... - [1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,8], ... - [-1.71, 1.71, 0.43, -8.75, 8.32, 8.32, -10.03, 1.71, 0.43, -1.12, 0.69, 0.035, -1.745, 1.71, ... - 0.43, -280 * y8 - 0.43, 280 * y8, -280 * y8, 0.43, 0.69, -1.81, 1.81, -280 * y6, 280 * y6, ... - -280 * y6]); - -% J2 = [ -% -1.71, 0.43, 8.32, 0, 0, 0, 0, 0; -% 1.71, -8.75, 0, 0, 0, 0, 0, 0; -% 0, 0, -10.03, 0.43, 0.035, 0, 0, 0; -% 0, 8.32, 1.71, -1.12, 0, 0, 0, 0; -% 0, 0, 0, 0, -1.745, 0.43, 0.43, 0; -% 0, 0, 0, 0.69, 1.71, -280 * y8 - 0.43, 0.69, -280 * y6; -% 0, 0, 0, 0, 0, 280 * y8, -1.81, 280 * y6; -% 0, 0, 0, 0, 0, -280 * y8, 1.81, -280 * y6]; +jac = sparse( ... + [1, 2, 1, 2, 4, 1, 3, 4, 3, 4, 6, 3, 5, 6, 5, 6, 7, 8, 5, 6, 7, 8, 6, 7, 8], ... + [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8], ... + [-k1, k1, k2, -k2 - k3, k3, k6, -k1 - k6, k1, k2, -k2 - k4, k4, k5, -k1 - k5, k1, k2, -k2 - kPlusE, kPlusE, ... + -kPlusE, k2, kMinus, -kSum, kSum, -kPlusPfrXPrime, kPlusPfrXPrime, -kPlusPfrXPrime]); end diff --git a/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m b/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m index 5f0260ba..cf80ea46 100644 --- a/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m +++ b/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m @@ -1,25 +1,17 @@ -function Jv = jacobianAdjointVectorProduct(~, y, v) +function vj = jacobianAdjointVectorProduct(~, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~) -y6 = conj(y(6)); -y8 = conj(y(8)); +kPlusPfrXPrime = conj(kPlus * y(6)); +kPlusE = conj(kPlus * y(8)); +kSum = conj(k2 + kMinus + kStar); -v1 = v(1); -v2 = v(2); -v3 = v(3); -v4 = v(4); -v5 = v(5); -v6 = v(6); -v7 = v(7); -v8 = v(8); - -Jv = [ - 1.71 * (v2 - v1); - 0.43 * v1 - 8.75 * v2 + 8.32 * v4; - 8.32 * v1 - 10.03 * v3 + 1.71 * v4; - 0.43 * v3 - 1.12 * v4 + 0.69 * v6; - 0.035 * v3 - 1.745 * v5 + 1.71 * v6; - 0.43 * (v5 - v6) + 280 * y8 * (v7 - v6 - v8); - 0.43 * v5 + 0.69 * v6 + 1.81 * (v8 - v7); - 280 * y6 * (v7 - v6 - v8)]; +vj = [ ... + conj(k1) * (v(2) - v(1)); + conj(k2) * v(1) - conj(k2 + k3) * v(2) + conj(k3) * v(4); + conj(k6) * v(1) - conj(k1 + k6) * v(3) + conj(k1) * v(4); + conj(k2) * v(3) - conj(k2 + k4) * v(4) + conj(k4) * v(6); + conj(k5) * v(3) - conj(k1 + k5) * v(5) + conj(k1) * v(6); + conj(k2) * v(5) - (conj(k2) + kPlusE) * v(6) + kPlusE * (v(7) - v(8)); + conj(k2) * v(5) + conj(kMinus) * v(6) + kSum * (v(8) - v(7)); + kPlusPfrXPrime * (v(7) - v(6) - v(8))]; end diff --git a/toolbox/+otp/+hires/jacobianVectorProduct.m b/toolbox/+otp/+hires/jacobianVectorProduct.m index df6095b3..a0f2c4f6 100644 --- a/toolbox/+otp/+hires/jacobianVectorProduct.m +++ b/toolbox/+otp/+hires/jacobianVectorProduct.m @@ -1,25 +1,28 @@ -function Jv = jacobianVectorProduct(~, y, v) +function jv = jacobianVectorProduct(~, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~) -y6 = y(6); -y8 = y(8); +k1Pr = k1 * v(1, :); +k2Pfr = k2 * v(2, :); +k3Pfr = k3 * v(2, :); +k1PrX = k1 * v(3, :); +k6PrX = k6 * v(3, :); +k2PfrX = k2 * v(4, :); +k4PfrX = k4 * v(4, :); +k1PrXPrime = k1 * v(5, :); +k5PrXPrime = k5 * v(5, :); +k2PfrXPrime = k2 * v(6, :); +k2PfrXPrimeE = k2 * v(7, :); +kMinusPfrXPrimeE = kMinus * v(7, :); +kPlusPfrXPrimeE = kPlus * (y(6) * v(8, :) + y(8) * v(6, :)); +jv7 = -k2PfrXPrimeE - (kMinus + kStar) * v(7, :) + kPlusPfrXPrimeE; -v1 = v(1); -v2 = v(2); -v3 = v(3); -v4 = v(4); -v5 = v(5); -v6 = v(6); -v7 = v(7); -v8 = v(8); - -Jv = [ - -1.71 * v1 + 0.43 * v2 + 8.32 * v3; - 1.71 * v1 - 8.75 * v2; - -10.03 * v3 + 0.43 * v4 + 0.035 * v5; - 8.32 * v2 + 1.71 * v3 - 1.12 * v4; - -1.745 * v5 + 0.43 * (v6 + v7); - 0.69 * (v4 + v7) + 1.71 * v5 - (280 * y8 + 0.43) * v6 - 280 * y6 * v8; - 280 * (y8 * v6 + y6 * v8) - 1.81 * v7; - 1.81 * v7 - 280 * (y8 * v6 + y6 * v8)]; +jv = [ ... + -k1Pr + k2Pfr + k6PrX; + k1Pr - k2Pfr - k3Pfr; + -k1PrX - k6PrX + k2PfrX + k5PrXPrime; + k3Pfr + k1PrX - k2PfrX - k4PfrX; + -k1PrXPrime - k5PrXPrime + k2PfrXPrime + k2PfrXPrimeE; + k4PfrX + k1PrXPrime - k2PfrXPrime + kMinusPfrXPrimeE - kPlusPfrXPrimeE; + jv7; + -jv7]; end From 94811f547970ad3d2a821e7b6bbad13b9f63aa80 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Sun, 13 Oct 2024 23:15:13 -0700 Subject: [PATCH 2/2] Finish HIRES docs --- docs/references.bib | 36 +++++++- toolbox/+otp/+hires/+presets/Canonical.m | 31 +++++++ toolbox/+otp/+hires/HIRESParameters.m | 20 ++--- toolbox/+otp/+hires/HIRESProblem.m | 89 ++++++++++++++++++- toolbox/+otp/+hires/f.m | 10 --- .../+hires/jacobianAdjointVectorProduct.m | 16 ++-- 6 files changed, 169 insertions(+), 33 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 728b7b8a..15a4c373 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -173,9 +173,39 @@ @article{Pet86 number = {4}, pages = {837--852}, publisher = {Society for Industrial and Applied Mathematics}, - title = {Order Results for Implicit Runge-Kutta Methods Applied to Differential/ Algebraic Systems}, + title = {Order Results for Implicit {Runge--Kutta} Methods Applied to Differential/Algebraic Systems}, volume = {23}, year = {1986}, - url = {https://www.jstor.org/stable/2157625} + doi = {10.1137/0723054} +} + +@article{Sch75, + title = {A new approach to explain the ``high irradiance responses'' of photomorphogenesis on the basis of phytochrome}, + volume = {2}, + ISSN = {1432-1416}, + DOI = {10.1007/bf00276015}, + number = {1}, + journal = {Journal of Mathematical Biology}, + author = {Schäfer, E.}, + year = {1975}, + pages = {41-56} +} + +@article{Got77, + title = {{MISS} - ein einfaches Simulations-System für biologische und chemische Prozesse}, + volume = {3}, + journal = {EDV in Medizin und Biologie}, + author = {B. A. Gottwald}, + year = {1977}, + pages = {85-90} +} + +@article{dSL98, + title={Collecting real-life problems to test solvers for implicit differential equations}, + author={de Swart, Jacques J. B. and Lioen, Walter M.}, + journal={CWI Quarterly}, + volume={11}, + number={1}, + pages={83-100}, + year={1998} } - diff --git a/toolbox/+otp/+hires/+presets/Canonical.m b/toolbox/+otp/+hires/+presets/Canonical.m index c4d1b49a..46f7b2c7 100644 --- a/toolbox/+otp/+hires/+presets/Canonical.m +++ b/toolbox/+otp/+hires/+presets/Canonical.m @@ -1,6 +1,37 @@ classdef Canonical < otp.hires.HIRESProblem + % HIRES preset proposed in :cite:p:`HW96` (pp. 144-145) which uses time span $t ∈ [0, 321.8122]$, initial conditions + % + % $$ + % y(0) = [1, 0, 0, 0, 0, 0, 0, 0.0057]^T, + % $$ + % + % and parameters + % + % $$ + % k_1 &= 1.71, \quad & k_2 &= 0.43, \quad & k_3 &= 8.32, \quad & k_4 &= 0.69, \quad & k_5 &= 0.035, \\ + % k_6 &= 8.32, & k_{+} &= 280, & k_{-} &= 0.69, & k^* &= 0.69, & o_{k_s} &= 0.0007. + % $$ + methods function obj = Canonical(varargin) + % Create the canonical HIRES problem object. + % + % Parameters + % ---------- + % varargin + % A variable number of name-value pairs. The accepted names are + % + % - ``K1`` – Value of $k_1$. + % - ``K2`` – Value of $k_2$. + % - ``K3`` – Value of $k_3$. + % - ``K4`` – Value of $k_4$. + % - ``K5`` – Value of $k_5$. + % - ``K6`` – Value of $k_6$. + % - ``KPlus`` – Value of $k_{+}$. + % - ``KMinus`` – Value of $k_{-}$. + % - ``KStar`` – Value of $k^*$. + % - ``OKS`` – Value of $o_{k_s}$. + tspan = [0, 321.8122]; y0 = [1; 0; 0; 0; 0; 0; 0; 0.0057]; params = otp.hires.HIRESParameters('K1', 1.71, 'K2', 0.43, 'K3', 8.32, 'K4', 0.69, 'K5', 0.035, ... diff --git a/toolbox/+otp/+hires/HIRESParameters.m b/toolbox/+otp/+hires/HIRESParameters.m index 5a8de173..82157a65 100644 --- a/toolbox/+otp/+hires/HIRESParameters.m +++ b/toolbox/+otp/+hires/HIRESParameters.m @@ -2,34 +2,34 @@ % Parameters for the HIRES problem. properties - % The stiffness parameter $ε$ for the "cusp catastrophe" model. + % The reaction rate $k_1$. K1 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_2$. K2 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_3$. K3 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_4$. K4 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_5$. K5 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_6$. K6 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_{+}$. KPlus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k_{-}$. KMinus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The reaction rate $k^*$. KStar %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} - % The diffusion coefficient $σ$ for all three variables. + % The source term $o_{k_s}$. OKS %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} end diff --git a/toolbox/+otp/+hires/HIRESProblem.m b/toolbox/+otp/+hires/HIRESProblem.m index 194592db..262ceb11 100644 --- a/toolbox/+otp/+hires/HIRESProblem.m +++ b/toolbox/+otp/+hires/HIRESProblem.m @@ -1,8 +1,71 @@ - classdef HIRESProblem < otp.Problem + % An eight-variable high irradiance response model from plant physiology. + % + % In :cite:p:`Sch75`, the following chemical reactions were proposed to explain the high irradiance responses of + % photomorphogenesis on the basis of phytochrome: + % + % $$ + % \begin{array}{cccc} + % \xrightarrow{o_{k_s}} & \ce{P_{{r}}} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}} \\ + % & {\scriptsize k_6} ↑ & & ↓ {\scriptsize k_3} \\ + % & \ce{P_{{r}}X} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X} \\ + % & {\scriptsize k_5} ↑ & & ↓ {\scriptsize k_4} \\ + % & \ce{P_{{r}}X'} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X'} \\ + % \end{array} + % $$ + % + % $$ + % \begin{array}{ccccc} + % \ce{E + P_{{r}}X'} & \xleftarrow{k_2} & \ce{P_{fr}X'E} & \xrightleftharpoons[k_{-}]{k_{+}} & \ce{P_{fr}X' + E} \\ + % & & ↓ {\scriptsize k^*} \\ + % & & \ce{P_{fr}' + X' + E} + % \end{array} + % $$ + % + % Reactants $\ce{P_{{r}}}$ and $\ce{P_{fr}}$ are the red and far-red absorbing form of phytochrome, respectively. + % These can be bound by receptors $\ce{X}$ and $\ce{X'}$, partially influenced by enzyme $\ce{E}$. The system is + % modeled by the differential equations :cite:p:`Got77` :cite:p:`dSL98` + % + % $$ + % \frac{d}{dt} \ce{[P_{{r}}]} &= -k_1 \ce{[P_{{r}}]} + k_2 \ce{[P_{fr}]} + k_6 \ce{[P_{{r}}X]} + o_{k_s}, \\ + % \frac{d}{dt} \ce{[P_{fr}]} &= k_1 \ce{[P_{{r}}]} - (k_2 + k_3) \ce{[P_{fr}]}, \\ + % \frac{d}{dt} \ce{[P_{{r}}X]} &= -(k_1 + k_6) \ce{[P_{{r}}X]} + k_2 \ce{[P_{fr}X]} + k_5 \ce{[P_{fr}X']}, \\ + % \frac{d}{dt} \ce{[P_{fr}X]} &= k_3 \ce{[P_{fr}]} + k_1 \ce{[P_{{r}}X]} - (k_2 + k_4) \ce{[P_{fr}X]}, \\ + % \frac{d}{dt} \ce{[P_{{r}}X']} &= -(k_1 + k_5) \ce{[P_{{r}}X']} + k_2 \ce{[P_{fr}X']} + k_2 \ce{[P_{fr}X'E]}, \\ + % \frac{d}{dt} \ce{[P_{fr}X']} &= k_4 \ce{[P_{fr}X]} + k_1 \ce{[P_{{r}}X']} - k_2 \ce{[P_{fr}X']} + % + k_{-} \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\ + % \frac{d}{dt} \ce{[P_{fr}X'E]} &= -(k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} + k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\ + % \frac{d}{dt} \ce{[E]} &= (k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}. + % $$ + % + % Notes + % ----- + % +---------------------+--------------------------------------------------------------------+ + % | Type | ODE | + % +---------------------+--------------------------------------------------------------------+ + % | Number of Variables | 8 | + % +---------------------+--------------------------------------------------------------------+ + % | Stiff | typically, depending on $k_1, …, k_6$, $k_{+}$, $k_{-}$, and $k^*$ | + % +---------------------+--------------------------------------------------------------------+ + % + % Example + % ------- + % >>> problem = otp.hires.presets.Canonical; + % >>> sol = problem.solve('AbsTol', 1e-12); + % >>> problem.plot(sol); methods function obj = HIRESProblem(timeSpan, y0, parameters) + % Create a HIRES problem object. + % + % Parameters + % ---------- + % timeSpan : numeric(1, 2) + % The start and final time. + % y0 : numeric(2, 1) + % The initial conditions. + % parameters : otp.hires.HIRESParameters + % The parameters. obj@otp.Problem('HIRES', 8, timeSpan, y0, parameters); end end @@ -26,7 +89,29 @@ function onSettingsChanged(obj) otp.hires.jacobianVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... 'JacobianAdjointVectorProduct', @(t, y, v) ... otp.hires.jacobianAdjointVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ... - 'Vectorized', 'on'); + 'Vectorized', 'on', ... + 'NonNegative', 1:obj.NumVars); + end + + function label = internalIndex2label(obj, index) + switch index + case 1 + label = 'P_r'; + case 2 + label = 'P_{fr}'; + case 3 + label = 'P_r X'; + case 4 + label = 'P_{fr} X'; + case 5 + label = 'P_r X'''; + case 6 + label = 'P_{fr} X'''; + case 7 + label = 'P_{fr} X'' E'; + case 8 + label = 'E'; + end end function fig = internalPlot(obj, t, y, varargin) diff --git a/toolbox/+otp/+hires/f.m b/toolbox/+otp/+hires/f.m index c7a275e0..691c5701 100644 --- a/toolbox/+otp/+hires/f.m +++ b/toolbox/+otp/+hires/f.m @@ -25,14 +25,4 @@ yPrime7; -yPrime7]; -% yPrime = [ -% -1.71 * y(1) + 0.43 * y(2) + 8.32 * y(3) + 0.0007; -% 1.71 * y(1) - 8.75 * y(2); -% -10.03 * y(3) + 0.43 * y(4) + 0.035 * y(5); -% 8.32 * y(2) + 1.71 * y(3) - 1.12 * y(4); -% -1.745 * y(5) + 0.43 * y(6) + 0.43 * y(7); -% -280 * y(6) * y(8) + 0.69 * y(4) + 1.71 * y(5) - 0.43 * y(6) + 0.69 * y(7); -% 280 * y(6) * y(8) - 1.81 * y(7); -% -280 * y(6) * y(8) + 1.81 * y(7)]; - end diff --git a/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m b/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m index cf80ea46..ad181b21 100644 --- a/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m +++ b/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m @@ -5,13 +5,13 @@ kSum = conj(k2 + kMinus + kStar); vj = [ ... - conj(k1) * (v(2) - v(1)); - conj(k2) * v(1) - conj(k2 + k3) * v(2) + conj(k3) * v(4); - conj(k6) * v(1) - conj(k1 + k6) * v(3) + conj(k1) * v(4); - conj(k2) * v(3) - conj(k2 + k4) * v(4) + conj(k4) * v(6); - conj(k5) * v(3) - conj(k1 + k5) * v(5) + conj(k1) * v(6); - conj(k2) * v(5) - (conj(k2) + kPlusE) * v(6) + kPlusE * (v(7) - v(8)); - conj(k2) * v(5) + conj(kMinus) * v(6) + kSum * (v(8) - v(7)); - kPlusPfrXPrime * (v(7) - v(6) - v(8))]; + conj(k1) * (v(2, :) - v(1, :)); + conj(k2) * v(1, :) - conj(k2 + k3) * v(2, :) + conj(k3) * v(4, :); + conj(k6) * v(1, :) - conj(k1 + k6) * v(3, :) + conj(k1) * v(4, :); + conj(k2) * v(3, :) - conj(k2 + k4) * v(4, :) + conj(k4) * v(6, :); + conj(k5) * v(3, :) - conj(k1 + k5) * v(5, :) + conj(k1) * v(6, :); + conj(k2) * v(5, :) - (conj(k2) + kPlusE) * v(6, :) + kPlusE * (v(7, :) - v(8, :)); + conj(k2) * v(5, :) + conj(kMinus) * v(6, :) + kSum * (v(8, :) - v(7, :)); + kPlusPfrXPrime * (v(7, :) - v(6, :) - v(8, :))]; end