Integration schemes for biochemical systems unconditional positivity and
Integration schemes for biochemical systems unconditional positivity and mass conservation Jorn Bruggeman Hans Burchard, Bob Kooi, Ben Sommeijer Theoretical Biology Vrije Universiteit, Amsterdam
Background Master Theoretical biology (2003) Start Ph. D study (2004) “Understanding the ‘organic carbon pump’ in mesoscale ocean flows” Focus: 1 D discretized water column turbulence and biota, simulation in time Tool: General Ocean Turbulence Model (GOTM) Modeling framework, split integration of advection, diffusion, production/destruction
Outline l Biochemical systems – – – l Traditional integration schemes – – l reaction-based framework conservation (of elements) positivity Euler, Runge-Kutta Modified Patankar New 1 st and 2 nd order schemes
Biochemical systems: the reaction • • chemical compounds = state variables c sources (left) are destroyed to produce sinks (right) constant stoichiometric coefficients (unit: compound/reaction) variable reaction rate (unit: reactions/time) Corresponding system of ODEs: Generalized for I state variables:
Systems of reactions Corresponding system of ODEs: Generalized for I state variables, R reactions :
The conservative reaction Compounds consist of chemical elements: Conservation: in reaction, no elements are created or destroyed! C O H for 1 conservative reaction:
Conservative systems With biochemical framework: microscopic conservation: in any reaction, no elements are created or destroyed Without biochemical framework: macroscopic conservation: in (closed) system, no elements are created or destroyed
Conservative integration schemes Macroscopic conservation: within system, quantities of element species are constant: Microscopic conservation? View on reaction-level is gone… ‘Biochemical integrity’: state variables change through known reactions only: for some vector If satisfied, implies microscopic/macroscopic conservation
Criteria for integration schemes Given a positive definite, conservative biochemical system: Integration scheme must satisfy: biochemical integrity/conservation: positivity: if given
Forward Euler, Runge-Kutta l l l Conservative: Non-positive Order: 1, 2, 4 etc.
Backward Euler, Gear l l Conservative: Positive for order 1 (Hundsdorfer & Verwer) Generalization to higher order eliminates positivity Slow! – – requires numerical approximation of partial derivatives requires solving linear system of equations
Modified Patankar: concepts l Burchard, Deleersnijder, Meister (2003) – l “A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations” Approach – – – Compound fluxes in production, destruction matrices (P, D) Pij = rate of conversion from j to i Dij = rate of conversion from i to j Source fluxes in D, sink fluxes in P
Modified Patankar: structure l l l Flux-specific multiplication factors cn+1/cn Represent ratio: (source after) : (source before) Multiple sources in reaction: – l multiple, different cn+1/cn factors Then: stoichiometric ratios not preserved!
Modified Patankar: example/conclusion l Conservative only if 1. 2. l l l every reaction contains ≤ 1 source compound source change ratios are identical (and remain so during simulation) Positive Order 1, 2 (higher possible? ) Requires solving linear system of equations
Typical MP conservation error Total nitrogen over 20 years: MP 1 st order MP-RK 2 nd order
New 1 st order scheme: structure l l Non-linear system of equations Positivity requirement fixes domain of product term p:
New 1 st order scheme: solution Component-wise, dividing by cn: Left and right, product over set Jn: l Polynomial for p: – l Derivative of polynomial < 0 within p domain: – l positive at left bound p=0, negative at right bound only one valid p Bisection technique is guaranteed to find p
New 1 st order scheme: conclusion l l l Positive Conservative: ± 20 bisection iterations (evaluations of polynomial) – – l Always cheaper than Backward Euler >4 state variables? Then cheaper than Modified Patankar Note: not suitable for stiff systems (unlike Modified Patankar)
Extension to 2 nd order
Test cases Linear system: Non-linear system:
Test case: linear system
Test case: non-linear system
Order tests Linear system: Non-linear system:
Plans l Publish new schemes – l Short term – – – l Bruggeman, Burchard, Kooi, Sommeijer (submitted 2005) Modeling ecosystems Aggregation into functional groups Modeling coagulation (marine snow) Extension to 3 D global circulation models
The end
- Slides: 25