Center Manifold and Lie Symmetry Calculations on a Quasi-chemical Model for Growth-death Kinetics in Food

TL;DR

This paper applies dynamical systems and Lie symmetry methods to analyze a bacterial growth-death model in food, providing insights into its mathematical structure and validating results against existing food science models.

Contribution

It introduces a novel application of center manifold and Lie symmetry calculations to a food microbiology model, enhancing understanding of its dynamics.

Findings

Calculation of center manifolds for the model

Identification of Lie symmetries under various parameters

Comparison with existing food science results

Abstract

Food scientists at the U.S. Army's Natick Solider Center have developed a model for the lifecyle of the bacteria Staphylococcus aureus in intermediate moisture bread. In this article, we study this model using dynamical systems and Lie symmetry methods. We calculate center manifolds and Lie symmetries for different cases of parameter values and compare our results to those of the food scientists.

Equations

$$v_1$$

$$v_2$$

$$v_3$$

$$v_4$$

$$\dot{M}$$

$$\dot{M^{\ast }}$$

$$\dot{A}$$

$$\dot{D}$$

$$\displaystyle\left(\begin{array}[]{c}y_{1}^{\prime}\\ y_{2}^{\prime}\\ y_{3}^{\prime}\end{array}\right)$$

$$\dot{x}$$

$$\dot{y}$$

$$\displaystyle\left(\begin{array}[]{c}y_{1}\\ y_{2}\\ y_{3}\end{array}\right)$$

$$u^{\prime }$$

$$\displaystyle\left(\begin{array}[]{c}v^{\prime}\\ w^{\prime}\end{array}\right)$$

$$u^{\prime }=0.$$

$$\displaystyle\left(\begin{array}[]{c}y_{1}^{\prime}\\ y_{2}^{\prime}\\ y_{3}^{\prime}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}y_{1}\\ y_{2}\\ y_{3}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}u^{\prime}\\ v^{\prime}\end{array}\right)$$

$$w^{\prime }$$

$$\dot{x}$$

$$\dot{y}$$

$$(Mh)(x)=Dh(x)[Ax+g(x,h(x))]-Bh(x)-j(x,h(x))=\overline{0}.$$

$$(Mh)(x)=\left(\begin{array}[]{c}-Gax^{2}+(-Gb-\varepsilon a+k_{1}\varepsilon c)x^{3}+O(x^{4})\\ k_{1}cx^{2}+k_{1}dx^{3}+O(x^{4})\end{array}\right).$$

$$(Mh)(x)$$

$$y_k^{\prime }={\omega }_k(t,y_1,y_2,\dots ,y_n),k=1,\dots ,n$$

$$X=\xi {\partial }_t+{\eta }_1{\partial }_{y_1}+{\eta }_2{\partial }_{y_2}+\cdots +{\eta }_n{\partial }_{y_n}$$

$$X^{(1)}(y_k^{\prime }-{\omega }_k)=0,k=1,\dots ,n$$

$$X^{(1)}=X+{\eta }_1^{(1)}{\partial }_{y_1^{\prime }}+{\eta }_2^{(1)}{\partial }_{y_2^{\prime }}+\cdots +{\eta }_n^{(1)}{\partial }_{y_n^{\prime }}$$

$${\eta }_k^{(1)}={\partial }_t{\eta }_k+y_1^{\prime }{\partial }_{y_1}{\eta }_k+y_2^{\prime }{\partial }_{y_2}{\eta }_k+\cdots +y_n^{\prime }{\partial }_{y_n}{\eta }_k-y_k^{\prime }({\partial }_t\xi +y_1^{\prime }{\partial }_{y_1}\xi +y_2^{\prime }{\partial }_{y_2}\xi +\cdots +y_n^{\prime }{\partial }_{y_n}\xi ).$$

$$X_1$$

$$X_2$$

$$X_3$$

$$X_4$$

$$X_5$$

$$X_6$$

$${\tilde{X}}_1$$

$${\tilde{X}}_2$$

$${\tilde{X}}_3$$

$${\tilde{X}}_4$$

Full text

1 Introduction

1.1 The model

To ensure the food that U.S. soldiers receive is as safe as possible, the growth of bacteria such as Staphylococcus aureus (S. aureus) needs to be addressed. The system of equations considered in this paper arises from a “quasi-chemical” kinetics model for the phases of the microbial life cycle of S. aureus in intermediate moisture bread. The food scientists who developed the model confirmed its usefulness by fitting to it the data from observations on bread crumbs with varying conditions of water activity, pH and temperature. The model differs from previous models in its attempt to model continuously the growth and death of the microorganism rather than focusing solely on either growth or inactivation. The model was developed by food scientists at the Natick Soldier Center, see Taub, et al  [4] for more on their techniques.

The model arose from the observation of four phases in the life cycle of S. aureus. The cells pass through the various stages of metabolizing ($M$), multiplying ($M^{*}$), sensitization to death ($M^{**}$), and dead ($D$). Additionally, the scientists hypothesized that there was an antagonist ($A$) present that would affect the cells. They found that without this added element their original model did not fit the observed data with any accuracy. The first step in the process describes cells moving from lag phase to growth phase ($M\rightarrow M^{*}$). In the next step, cells multiply via binary division and then the newly multiplied cells interact with an antagonist ($M^{*}\rightarrow 2M^{*}+A$). The last two steps represent two different pathways to death: the first with cells interacting with an antagonist, then passing to sensitization before death ($A+M^{*}\rightarrow M^{**}\rightarrow D)$ and lastly the cells experiencing natural death ($M^{*}\rightarrow D$).

The following equations represent the velocities of each of the above steps ($v$) as they relate to the concentrations of cells in various the phases. Each equation has a rate constant ($k$) associated to it.

[TABLE]

Finally these velocities are represented by the following system of ordinary differential equations:

[TABLE]

where $G=k_{2}-k_{4}$ is the net natural growth rate and $\varepsilon=10^{-9}k_{3}$. It is assumed that all the rate constants have non-negative values. The initial conditions at time zero are $M(0)=I$, the inoculum level $I\approx 10^{3}-10^{4}$, and $M^{*}(0)=A(0)=D(0)=0$.

1.2 A simplification

We notice that the fourth equation is uncoupled since there are no terms involving the variable $D$ in any of the other equations and $\dot{D}$ depends on $M^{*}$ and $A$. Therefore to investigate the dynamics of our system, we reduce to a system of three equations. Renaming our variables ($y_{1}=M,\ y_{2}=M^{*},\ y_{3}=A$) we have the following system equivalent to equations 5-8:

[TABLE]

2 Normalizing the system

To consider the invariant manifold structure of a system, it is necessary to write the system in normal form as follows:

[TABLE]

with $(x,y)\in\mathbb{R}^{n}\times\mathbb{R}^{m}$, the $n\times n$ matrix $A$ having eigenvalues with zero real part and the $m\times m$ matrix $B$ having eigenvalues with nonzero real part. The functions $g(x,y)$ and $j(x,y$) must be zero with zero first partial derivatives at the origin.

The system (21) above is not in normal form since the $y_{3}^{\prime}$ equation corresponds to the zero eigenvalue piece and the nonlinear term of $y_{3}^{\prime}$ does not have all zero partial derivatives at the origin. Thus we must normalize by a change of coordinates using the eigenvectors of the matrix of the linear terms of the equation. We will investigate the invariant manifolds in a neighborhood of $G=0$. Writing our system in normal form for nonzero $G$ does not depend on the sign of $G$, so we treat the negative and positive case simultaneously. We let $T$ be the matrix of eigenvectors of the eigenvalues of the linear terms of our system and let

[TABLE]

Using the inverse of the matrix $T$ we can solve for $u,\ v$ and $w$, find their derivatives and finally write our system in normal form as follows:

[TABLE]

where $f(u,v,w)=-\frac{\varepsilon}{G}(v-k_{1}w)(k_{2}v+G(u+k_{2}w))$. Since $f(u,v,w)$ and its first partials with respect to $u,\ v,$ and $w$ are all zero at the origin $(u,v,w)=(0,0,0)$, we have our system in normal form and we see immediately that we have a one-dimensional center manifold in the case that $G\neq 0$. For $G>0$, we also have a one-dimensional stable and a one-dimensional unstable manifold. For $G<0$, we have a two-dimensional stable manifold. The system reduced to the center manifold simply becomes

[TABLE]

For the case $G=0$, we have a slightly simpler system of equations:

[TABLE]

Note that we now have two zero eigenvalues and one negative eigenvalue for the matrix in the linear term. Since zero is a repeated eigenvalue, we must use generalized eigenvectors to find the normalization of this system. Three such eigenvectors are $(0,0,1),\ (0,1,0),$ and $(1,-1,\frac{k_{2}}{k_{1}})$. Then to transform our system we again let $T$ be the matrix consisting of these eigenvectors and let

[TABLE]

As above, this allows us to write our system in normal form:

[TABLE]

where $g(u,v,w)=-\varepsilon(v-w)(u+\frac{k_{2}}{k_{1}}w)$. Since $g(u,v,w)$ is zero at the origin and all of its first partial derivatives are also zero at the origin, we can see that we have a two dimensional center manifold and a one dimensional stable manifold.

3 Center manifold calculations

Recall that a center manifold $W^{c}=\{(x,y)|y=h(x)\}$ is described by $h(x)$ where $h(0)=Dh(0)=0$. We consider a system written in normal form

[TABLE]

with $A$ having eigenvalues with zero real part and $B$ eigenvalues with nonzero real part. Then we determine $h(x)$ by finding the function that satisfies the following condition:

[TABLE]

The sign of $G$ does not change the outcome of this calculation, thus we treat the case $G\neq 0$ at once. We have $h:V\rightarrow\mathbb{R}^{2}$, $V\subset\mathbb{R}$ a neighborhood of the origin. Thus let $h(x)=(h_{1}(x),h_{2}(x))=(ax^{2}+bx^{3}+O(x^{4}),cx^{2}+dx^{3}+O(x^{4}))$. Then $f(x,h_{1}(x),h_{2}(x))=\varepsilon(-a+k_{1}c)x^{3}+O(x^{4})$ resulting in

[TABLE]

Solving for $(Mh)(x)=\overline{0}$, $h_{1}(x)=h_{2}(x)=O(x^{4})$. Thus up to third order, we have $h_{1}(x)=h_{2}(x)=0$, so a center manifold is simply the $u-$axis.

Next we consider the case $G=0$. Here $h:V\rightarrow\mathbb{R}$, $V\subset\mathbb{R}^{2}$, a neighborhood of the origin. We let $h(x)=h(x_{1},x_{2})=ax_{1}^{2}+bx_{2}^{2}+cx_{1}x_{2}+dx_{1}^{3}+ex_{2}^{3}+fx_{1}^{2}x_{2}+jx_{1}x_{2}^{2}$. Then we calculate

[TABLE]

resulting in $h(x_{1},x_{2})=O(x^{4})$, thus $h(x_{1},x_{2})=0$ up to order three. Hence in this case the $uv-$plane is a center manifold.

4 Lie Symmetry

Recall that a Lie symmetry is a map from the set of solutions of a system of differential equations to the set itself. For a system of first order ordinary differential equations

[TABLE]

the Lie symmetries that transform the variables $t,y_{1},\dots,y_{n}$ have infinitesimal generators of the form

[TABLE]

where $\xi=\xi(t,y_{1},y_{2},\dots,y_{n})$ and $\eta_{k}=\eta_{k}(t,y_{1},y_{2},\dots,y_{n})$ for all $k$. The infinitesimal generator must satisfy the Linearized Symmetry Condition:

[TABLE]

when (87) holds. In this case the prolongation of $X$ is as follows:

[TABLE]

where $\eta_{k}^{(1)}$ is defined as $\eta_{k}^{(1)}=D_{t}\eta_{k}-y_{k}^{\prime}D_{t}\xi$. The total derivative $D_{t}$ in this case is $D_{t}=\partial_{t}+y_{1}^{\prime}\partial_{y_{1}}+\cdots+y_{n}^{\prime}\partial_{y_{n}}$. Thus we have the following:

[TABLE]

A system of first order ODEs has an infinite number of symmetries. We find symmetries by solving for the functions $\xi,\eta_{k}$ that satisfy the Linearized Symmetry Condition (89). This condition reduces to a system of PDEs which are computationally difficult to solve. We use the “Intro to Symmetry” package in Mathematica and a script included in Cantwell [1] to calculate the symmetries for our system. We are limited in the symmetries we can calculate by our computing power. In the case $G\neq 0$ we calculate symmetries up to third order in our original coordinates $y_{1},\ y_{2},$ and $y_{3}$ and then use a change of coordinates on our symmetries to rewrite in the coordinates $u,\ v,$ and $w$ of our equations in normal form. Since the case $G=0$ involves simpler equations, we are able to calculate these symmetries directly from the equations in normal form, however we followed the same method as in the $G\neq 0$ case since we want to be able to compare cases.

4.1 The case $G\neq 0$

The infinitesimals of the Lie symmetries (up to order 3) are listed in an array with $\{\xi,\eta_{1},\eta_{2},\eta_{3}\}$, representing the infinitesimal generator $X=\xi\partial_{t}+\eta_{1}\partial_{y_{1}}+\eta_{2}\partial_{y_{2}}+\eta_{3}\partial_{y_{3}}$.

[TABLE]

Then we transform the infinitesimal generators of the Lie symmetries found in the $y_{i}$ coordinates as follows. If $X$ is an infinitesimal generator in $y_{i}$, then $\tilde{X}=(Xt)\partial_{t}+(Xu)\partial_{u}+(Xv)\partial_{v}+(Xw)\partial_{w}$ is the corresponding infinitesimal generator for a Lie symmetry in the $u,v,w$ coordinates The transformed symmetries in the form $\tilde{X}=\{\tilde{\xi},\tilde{\eta_{1}},\tilde{\eta_{2}},\tilde{\eta_{3}}\}$ where $\tilde{X}=\tilde{\xi}\partial_{t}+\tilde{\eta_{1}}\partial_{u}+\tilde{\eta_{2}}\partial_{v}+\tilde{\eta_{3}}\partial_{w}$:

[TABLE]

where $f(u,v,w)$ is as above, $j(u,v,w)=v-k_{1}w$, $l(u,v,w)=u+\frac{k_{2}}{G}v+k_{2}w$ and $m(u,v,w)=-\frac{1}{\varepsilon}(G+k_{1})w$.

4.2 The case $G=0$

Again we calculate the infinitesimals of the Lie symmetries (up to order 3) of the original system with coordinates $\{y_{1},y_{2},y_{3}\}$ and list them as $X=\{\xi,\eta_{1},\eta_{2},\eta_{3}\}$, representing the infinitesimal generator $X=\xi\partial_{t}+\eta_{1}\partial_{y_{1}}+\eta_{2}\partial_{y_{2}}+\eta_{3}\partial_{y_{3}}$.

[TABLE]

Then we transform these to the $u,v,w$ coordinate system as above with $\tilde{X}=\{\tilde{\xi},\tilde{\eta_{1}},\tilde{\eta_{2}},\tilde{\eta_{3}}\}$ where $\tilde{X}=\tilde{\xi}\partial_{t}+\tilde{\eta_{1}}\partial_{u}+\tilde{\eta_{2}}\partial_{v}+\tilde{\eta_{3}}\partial_{w}$:

[TABLE]

where $n(u,v,w)=\varepsilon k_{2}w(-v+w)+k_{1}(k_{2}v+\varepsilon u(-v+w))$ and $p(u,v,w)=(v-w)(u+\frac{k_{2}}{k_{1}}w)$.

5 The connections between the center manifold and the Lie symmetry

Recently Cicogna and Gaeta [2] have written about the connections between dynamical systems and Lie symmetries. We are interested in particular in their results on invariant manifolds. They have commented that any Lie symmetry of the system will leave invariant both the stable and unstable manifolds. Due to the non-uniqueness of center manifolds, a Lie symmetry will map a center manifold to another (possibly the same) center manifold. The following result indicates when a center manifold given by $\omega(u)$ will be invariant under a given Lie symmetry, in their notation $X=\phi\partial_{u}+\psi\partial_{v}$.

Lemma 5.1** (Lemma 4 of [2] Chapter 7).**

A center manifold $w(u)$ is invariant under a Lie symmetry $X=\phi\partial_{u}+\psi\partial_{v}$ if and only if

[TABLE]

For the case $G\neq 0$, $\omega(u)=\{0,0\}$ giving zero on the right side of this equality. Thus the left side of this equation evaluated on the center manifold must always be zero if our center manifold is to be invariant under the action of the symmetry. This is the case with all of our Lie symmetries as given above. For example consider $X_{2}$ with $\phi(u,v,w)=\frac{1}{G}(G-k_{2})j(u,v,w)f(u,v,w)$ and $\psi(u,v,w)=\{j(u,v,w)(Gv+f(u,v,w)),-k_{1}wj(u,v,w)\}$. Since $j(u,0,0)\equiv 0$, $\psi(u,\omega(u))=\psi(u,0,0)=\{0,0\}$, thus satisfying the necessary and sufficient condition of the lemma. It is easy to determine that the remainder of the symmetries in this case also leave the center manifold invariant. Thus the center manifolds inherit these Lie symmetries. However, in this case, since $v=w=0$, all of our symmetries become trivial.

Recall that in the case $G=0$ we found a center manifold to be the $uv-$plane. Now, in the notation of our lemma, $\omega(u)=0$, and again the right side of our equation is zero. Thus we must have $\psi(u,v,0)=0$ for any symmetry that leaves invariant this center manifold. It can be easily checked to see that all of the symmetries listed above do indeed satisfy this necessary and sufficient condition. In this case the center manifold again inherits the Lie symmetries which are now nontrivial. The restriction of the system to our center manifold, the $uv-$plane, is

[TABLE]

The nontrivial symmetries inherited by this system are

[TABLE]

If we transform back to our original variables, we see that on the center manifold $u=y_{3}$ and $v=y_{2}$, resulting in the system:

[TABLE]

and the symmetries:

[TABLE]

While we have calculated the infinitesimal generators, it would be interesting to determine the actual Lie symmetries on the center manifolds. We would like to say precisely what these maps do to various trajectories on the center manifold and to the flow in general. This is however, a very difficult question. There is no known method that allows us to take the infinitesimal generators of any Lie symmetry and integrate them to find the actual symmetries. The difficulty of this question is analogous to the solving of a system of differential equations analytically.

For example, if we consider $\hat{X}_{4}$ with $\eta_{2}(t,y_{2},y_{3})=\frac{-k_{2}}{\varepsilon}y_{2}+y_{2}y_{3}$ and $\eta_{3}=y_{2}y_{3}$, this means that, letting $\gamma$ be the parameter of the one-parameter Lie group, we need to solve the following for $\hat{y_{2}}$ and $\hat{y_{3}}$, giving us the map $(\hat{y_{2}},\hat{y_{3}})$ as our symmetry:

[TABLE]

This is equivalent to the system above. Attempting to solve this system we find it equivalent to solving the following:

[TABLE]

with the initial conditions $\hat{y_{2}}(\gamma,y_{2},y_{3})|_{\gamma=0}=y_{2}$ and $\hat{y_{3}}(\gamma,y_{2},y_{3})|_{\gamma=0}=y_{3}$.

This is something we continue to work on for this particular system as well as in general.

5.1 Comparison to previous results

Based on numerical solutions of the original system of equations Ross et. al [3] predicted trajectories for $M,\ M^{*},A$ and $D$ with particular emphasis on the concentrations of $M^{*}$ (cells undergoing multiplication) and $A$ (the antagonist). They found that the behavior depended on the values of the various constants $k_{i}$. In particular, with $k_{3}=0$ and $G>0$, they found unrestrained growth of both $M^{*}$ and $A$. For the values $k_{3}=0$ and a negative $G$, $A$ increases toward an upper limit and $M^{*}$ increases slightly but then begins to decrease toward zero. For $k_{3}>0$ and $G>0$, both $M^{*}$ and $A$ increase initially, but then $M^{*}$ reaches a maximum and begins to decline while $A$ approaches an upper bound. All of these analyses combined to indicate to the food scientists that the necessary constraints for growth-death kinetics are non-zero values for $k_{3}$ and positive values of $G$.

In our consideration of the system, we also found that $M^{*}(=y_{3})$ and $A(=y_{2})$ were the two variables that determined the behavior of the system. In the $G\neq 0$ case, the center manifold is the $u$-axis, which corresponds to $A$ when all other variables are zero, as on the center manifold. When $G=0$, the reduced system on the center manifold is given by equations 92 and 93. An inspection of this system, noting that $u=A$ and $v=M^{*}$, shows that the behavior is qualitatively identical to that found numerically in [4] for the case $k=[1\ 4\ 100\ 4]$, i.e. $k_{2}=k_{4}=4$ resulting in $G=0$. In both the results are that $M^{*}$ goes to zero and $A$ approaches a constant value.

Acknowledgements

This research was performed while the first author held a National Research Council Research Associateship Award jointly at the U.S. Army Natick Soldier Center, Natick, Massachusetts and the United States Military Academy, West Point, New York.