Center Manifold and Lie Symmetry Calculations
on a Quasi-chemical Model for Growth-death
Kinetics in Food
Rachelle C. DeCoste
Department of Mathematical Sciences
The United States Military Academy
West Point, NY 10996
Tel.: 845-938-2530
Rachelle.DeCoste@usma.edu
Louis Piscitelle
U.S. Army RDECOM, Natick Soldier Center
Kansas Street, Natick, MA 01760
Tel.: 508-233-4294
Louis.Piscitelle@us.army.mil
November 28, 2018
Abstract
Food scientists at the U.S. Army’s Natick Solider Center have devel-
oped a model for the lifecyle of the bacteria Staphylococcus aureus in
intermediate moisture bread. In this article, we study this model us-
ing 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.
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 intermedi-
ate moisture bread. The food scientists who developed the model confirmed its
usefulness by fitting to it the data from observations on bread crumbs with vary-
ing conditions of water activity, pH and temperature. The model differs from
http://arxiv.org/abs/0704.1284v1
CENTER MANIFOLDS AND LIE SYMMETRIES 2
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), mul-
tiplying (M∗), sensitization to death (M∗∗), and dead (D). Additionally, the
scientists hypothesized that there was an antagonist (A) present that would af-
fect 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 →M∗). In the next
step, cells multiply via binary division and then the newly multiplied cells in-
teract with an antagonist (M∗ → 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∗ →M∗∗ → D) and lastly the
cells experiencing natural death (M∗ → 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.
v1 = k1M (1)
v2 = k2M
∗ (2)
v3 = (10
−9)k3M
∗A (3)
v4 = k4M
∗ (4)
Finally these velocities are represented by the following system of ordinary dif-
ferential equations:
Ṁ = −v1 = −k1M (5)
Ṁ∗ = v1 + v2 − v3 − v4 = k1M +M
∗(G− εA) (6)
Ȧ = v2 − v3 =M
∗(k2 − εA) (7)
Ḋ = v3 + v4 =M
∗(k4 + εA) (8)
where G = k2−k4 is the net natural growth rate and ε = 10
−9k3. It is assumed
that all the rate constants have non-negative values. The initial conditions at
time zero areM(0) = I, the inoculum level I ≈ 103− 104, 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 involv-
ing the variable D in any of the other equations and 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 (y1 =M, y2 = M
∗, y3 = A) we have
CENTER MANIFOLDS AND LIE SYMMETRIES 3
the following system equivalent to equations 5-8:
−k1 0 0
k1 G 0
0 k2 0
−εy2y3
−εy2y3
 (9)
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:
ẋ = Ax+ g(x, y) (10)
ẏ = By + j(x, y) (11)
with (x, y) ∈ Rn × Rm, the n × n matrix A having eigenvalues with zero real
part and the m×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 (9) above is not in normal form since the y′3 equation corresponds
to the zero eigenvalue piece and the nonlinear term of y′3 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
 = T
0 0 G+ k1
0 1 −k1
 . (12)
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:
u′ = 0 · u+
f(u, v, w) (13)
0 −k1
f(u, v, w)
where f(u, v, w) = − ε
(v−k1w)(k2v+G(u+k2w)). 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 6= 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
u′ = 0. (15)
CENTER MANIFOLDS AND LIE SYMMETRIES 4
For the case G = 0, we have a slightly simpler system of equations:
−k1 0 0
k1 0 0
0 k2 0
−εy2y3
−εy2y3
 (16)
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, k2
). Then to transform our
system we again let T be the matrix consisting of these eigenvectors and let
 = T
0 0 1
0 1 −1
1 0 k2
 . (17)
As above, this allows us to write our system in normal form:
g(u, v, w)
g(u, v, w)
w′ = −k1w + 0 (19)
where g(u, v, w) = −ε(v−w)(u+ k2
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 manifoldW 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
ẋ = Ax+ g(x, y) (20)
ẏ = By + j(x, y) (21)
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:
(Mh)(x) = Dh(x)[Ax + g(x, h(x))] −Bh(x)− j(x, h(x)) = 0.
The sign of G does not change the outcome of this calculation, thus we treat
the case G 6= 0 at once. We have h : V → R2, V ⊂ R a neighborhood of the
origin. Thus let h(x) = (h1(x), h2(x)) = (ax
2+ bx3+O(x4), cx2+dx3+O(x4)).
Then f(x, h1(x), h2(x)) = ε(−a+ k1c)x
3 +O(x4) resulting in
(Mh)(x) =
−Gax2 + (−Gb− εa+ k1εc)x
3 +O(x4)
2 + k1dx
3 +O(x4)
CENTER MANIFOLDS AND LIE SYMMETRIES 5
Solving for (Mh)(x) = 0, h1(x) = h2(x) = O(x
4). Thus up to third order, we
have h1(x) = h2(x) = 0, so a center manifold is simply the u−axis.
Next we consider the case G = 0. Here h : V → R, V ⊂ R2, a neighborhood
of the origin. We let h(x) = h(x1, x2) = ax
2+cx1x2+dx
2. Then we calculate
(Mh)(x) =
hx1(x1, x2), hx2(x1, x2)
k2x2 + g(x1, x2, h(x1, x2))
g(x1, x2, h(x1, x2))
+ k1h(x1, x2)
= (2ak2 + ck1)x1x2 + (ak1)x
1 + (ck2 + bk1)x
2 + (dk1)x
1 + (jk2 + ek1)x
+(3dk2 − 2aε− cε+ fk1)x
1x2 + (2fk2 − cε− 2bε+ jk1)x1x
resulting in h(x1, x2) = O(x
4), thus h(x1, x2) = 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
k = ωk(t, y1, y2, . . . , yn), k = 1, . . . , n (22)
the Lie symmetries that transform the variables t, y1, . . . , yn have infinitesimal
generators of the form
X = ξ∂t + η1∂y1 + η2∂y2 + · · ·+ ηn∂yn (23)
where ξ = ξ(t, y1, y2, . . . , yn) and ηk = ηk(t, y1, y2, . . . , yn) for all k. The in-
finitesimal generator must satisfy the Linearized Symmetry Condition:
(1)(y′k − ωk) = 0, k = 1, . . . , n (24)
when (22) holds. In this case the prolongation of X is as follows:
X(1) = X + η
1 ∂y′1 + η
2 ∂y′2 + · · ·+ η
n ∂y′n (25)
where η
is defined as η
= Dtηk − y
kDtξ. The total derivative Dt in this
case is Dt = ∂t + y
1∂y1 + · · ·+ y
n∂yn . Thus we have the following:
= ∂tηk+y
1∂y1ηk+y
2∂y2ηk+· · ·+y
n∂ynηk−y
k(∂tξ+y
1∂y1ξ+y
2∂y2ξ+· · ·+y
n∂ynξ).
A system of first order ODEs has an infinite number of symmetries. We
find symmetries by solving for the functions ξ, ηk that satisfy the Linearized
Symmetry Condition (24). 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
CENTER MANIFOLDS AND LIE SYMMETRIES 6
for our system. We are limited in the symmetries we can calculate by our
computing power. In the case G 6= 0 we calculate symmetries up to third order
in our original coordinates y1, y2, and y3 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 6= 0 case since we want to be able to
compare cases.
4.1 The case G 6= 0
The infinitesimals of the Lie symmetries (up to order 3) are listed in an array
with {ξ, η1, η2, η3}, representing the infinitesimal generator X = ξ∂t + η1∂y1 +
η2∂y2 + η3∂y3 .
X1 = {1, 0, 0, 0}
X2 = {y2,−k1y1y2, k1y1y2 +Gy
2 − εy
2y3, k2y
2 − εy
X3 = {y3,−k1y1y3, k1y1y3 +Gy2y3 − εy2y
3 , k2y2y3 − εy2y
X4 = {0,−y1, y1 +
y2y3,
y2y3}
X5 = {
t,−ty1, ty1 +
ty2 −
ty2y3,
ty2 −
ty2y3}
X6 = {
y1y2 + y1y2y3,
y1y2 + y1y2y3}
Then we transform the infinitesimal generators of the Lie symmetries found in
the yi coordinates as follows. If X is an infinitesimal generator in yi, then X̃ =
(Xt)∂t+(Xu)∂u+(Xv)∂v+(Xw)∂w is the corresponding infinitesimal generator
for a Lie symmetry in the u, v, w coordinates The transformed symmetries in
the form X̃ = {ξ̃, η̃1, η̃2, η̃3} where X̃ = ξ̃∂t + η̃1∂u + η̃2∂v + η̃3∂w:
X̃1 = {1, 0, 0, 0}
X̃2 = {j(u, v, w),
(G− k2)j(u, v, w)f(u, v, w), j(u, v, w)(Gv + f(u, v, w)),−k1wj(u, v, w)}
X̃3 = {l(u, v, w),
(G− k2)l(u, v, w)f(u, v, w), l(u, v, w)(Gv + f(u, v, w)),−k1wl(u, v, w)}
X̃4 = {0,
(G− k2)f(u, v, w),
(Gv + f(u, v, w)),−w}
X̃5 = {
(G− k2)f(u, v, w),
(Gv + f(u, v, w)),−tw}
X̃6 = {m(u, v, w),
(G− k2)m(u, v, w)f(u, v, w),m(u, v, w)(Gv + f(u, v, w)),−k1wm(u, v, w)
where f(u, v, w) is as above, j(u, v, w) = v − k1w, l(u, v, w) = u +
v + k2w
and m(u, v, w) = − 1
(G+ k1)w.
CENTER MANIFOLDS AND LIE SYMMETRIES 7
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 {y1, y2, y3} and list them as X = {ξ, η1, η2, η3},
representing the infinitesimal generator X = ξ∂t + η1∂y1 + η2∂y2 + η3∂y3 .
X1 = {1, 0, 0, 0}
X2 = {y2,−k1y1y2, k1y1y2 − εy
2y3, k2y
2 − εy
X3 = {y3,−k1y1y3, k1y1y3 − εy2y
3 , k2y2y3 − εy2y
X4 = {0,
y1 + y2y3,
y2 + y2y3}
X5 = {
ty1 + ty2y3,
ty2 + ty2y3}
X6 = {
1 + y1y2y3,
y1y2 + y1y2y3}
Then we transform these to the u, v, w coordinate system as above with X̃ =
{ξ̃, η̃1, η̃2, η̃3} where X̃ = ξ̃∂t + η̃1∂u + η̃2∂v + η̃3∂w:
X̃1 = {1, 0, 0, 0}
X̃2 = {v − w,
(v − w)n(u, v, w),−ε(v − w)p(u, v, w),−k1(v − w)w}
X̃3 = {u+
w)n(u, v, w),−ε(u +
w)p(u, v, w),−k1(u+
X̃4 = {0,
n(u, v, w), p(u, v, w),
X̃5 = {
tn(u, v, w), tp(u, v, w),
X̃6 = {
wn(u, v, w), wp(u, v, w),
where n(u, v, w) = εk2w(−v + w) + k1(k2v + εu(−v + w)) and p(u, v, w) =
(v − w)(u + k2
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 ω(u) will be invariant under a given
Lie symmetry, in their notation X = φ∂u + ψ∂v.
CENTER MANIFOLDS AND LIE SYMMETRIES 8
Lemma 5.1 (Lemma 4 of [2] Chapter 7). A center manifold w(u) is invariant
under a Lie symmetry X = φ∂u + ψ∂v if and only if
ψ(u, ω(u)) = (∂u(ω(u))) · φ(u, ω(u)).
For the case G 6= 0, ω(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 X2 with φ(u, v, w) =
(G − k2)j(u, v, w)f(u, v, w) and
ψ(u, v, w) = {j(u, v, w)(Gv + f(u, v, w)),−k1wj(u, v, w)}. Since j(u, 0, 0) ≡ 0,
ψ(u, ω(u)) = ψ(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, ω(u) = 0, and again the right side of our
equation is zero. Thus we must have ψ(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
u′ = k2v − εuv (27)
v′ = −εuv. (28)
The nontrivial symmetries inherited by this system are
X̂2 = {v,
vn(u, v, 0),−εvp(u, v, 0), 0}
X̂3 = {u,
un(u, v, 0),−εup(u, v, w), 0}
X̂4 = {0,
n(u, v, 0), p(u, v, 0), 0}
X̂5 = {
tn(u, v, 0), tp(u, v, 0), 0}
If we transform back to our original variables, we see that on the center
manifold u = y3 and v = y2, resulting in the system:
y′2 = −εy2y3 (29)
y′3 = k2y2 − εy2y3 (30)
CENTER MANIFOLDS AND LIE SYMMETRIES 9
and the symmetries:
X̂2 = y2∂t +
2 − εy
∂y2 − εy
2y3∂y3 (31)
X̂3 = y3∂t +
k2y2y3 − εy2y
∂y2 − εy2y
3∂y3 (32)
X̂4 =
y2 + y2y3
∂y2 + y2y3∂y3 (33)
X̂5 =
t∂t +
ty2 + ty2y3
∂y2 + ty2y3∂y3 (34)
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 X̂4 with η2(t, y2, y3) =
y2 + y2y3 and η3 =
y2y3, this means that, letting γ be the parameter of the one-parameter Lie
group, we need to solve the following for ŷ2 and ŷ3, giving us the map (ŷ2, ŷ3)
as our symmetry:
ŷ2 + ŷ2ŷ3 (35)
= ŷ2ŷ3 (36)
This is equivalent to the system above. Attempting to solve this system we find
it equivalent to solving the following:
ŷ2 = e
+ŷ3)dγ (37)
ŷ3 = e
ŷ2dγ (38)
with the initial conditions ŷ2(γ, y2, y3)|γ=0 = y2 and ŷ3(γ, y2, y3)|γ=0 = y3.
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
ki. In particular, with k3 = 0 and G > 0, they found unrestrained growth of
both M∗ and A. For the values k3 = 0 and a negative G, A increases toward an
upper limit and M∗ increases slightly but then begins to decrease toward zero.
CENTER MANIFOLDS AND LIE SYMMETRIES 10
For k3 > 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 k3 and positive
values of G.
In our consideration of the system, we also found thatM∗(= y3) and A(= y2)
were the two variables that determined the behavior of the system. In the G 6= 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 27 and 28. 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.
k2 = k4 = 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.
References
[1] B. J. Cantwell, Introduction to Symmetry Analysis, Cambridge Univer-
sity Press, Cambridge, United Kingdom, 2002.
[2] G. Cicogna and G. Gaeta, Symmetry and Perturbation Theory in Non-
linear Dynamics, Springer-Verlag, 1999.
[3] E. Ross, I. Taub, C. Doona, F. Feeherry, K. Kustin, The math-
ematical properties of the quasi-chemical model for microorganism growth
– death kinetics in food, International Journal of Food Microbiology, 99
(2005), pp. 157–171.
[4] I. A. Taub, F. E. Feeherry, E. W. Ross, K. Kustin, and
C. J.Doona, A Quasi-Chemical Kinetics Model for the Growth and Death
of Staphylococcus aureus in Intermediate Moisture Bread, Journal of Food
Science, 68, No. 8 (2003), pp. 2530–2537.
	Introduction
	The model
	A simplification
	Normalizing the system
	Center manifold calculations
	Lie Symmetry
	The case G=0
	The case G=0
	The connections between the center manifold and the Lie symmetry
	Comparison to previous results

  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.

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 intermedi-
ate moisture bread. The food scientists who developed the model confirmed its
usefulness by fitting to it the data from observations on bread crumbs with vary-
ing conditions of water activity, pH and temperature. The model differs from
http://arxiv.org/abs/0704.1284v1
CENTER MANIFOLDS AND LIE SYMMETRIES 2
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), mul-
tiplying (M∗), sensitization to death (M∗∗), and dead (D). Additionally, the
scientists hypothesized that there was an antagonist (A) present that would af-
fect 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 →M∗). In the next
step, cells multiply via binary division and then the newly multiplied cells in-
teract with an antagonist (M∗ → 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∗ →M∗∗ → D) and lastly the
cells experiencing natural death (M∗ → 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.
v1 = k1M (1)
v2 = k2M
∗ (2)
v3 = (10
−9)k3M
∗A (3)
v4 = k4M
∗ (4)
Finally these velocities are represented by the following system of ordinary dif-
ferential equations:
Ṁ = −v1 = −k1M (5)
Ṁ∗ = v1 + v2 − v3 − v4 = k1M +M
∗(G− εA) (6)
Ȧ = v2 − v3 =M
∗(k2 − εA) (7)
Ḋ = v3 + v4 =M
∗(k4 + εA) (8)
where G = k2−k4 is the net natural growth rate and ε = 10
−9k3. It is assumed
that all the rate constants have non-negative values. The initial conditions at
time zero areM(0) = I, the inoculum level I ≈ 103− 104, 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 involv-
ing the variable D in any of the other equations and 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 (y1 =M, y2 = M
∗, y3 = A) we have
CENTER MANIFOLDS AND LIE SYMMETRIES 3
the following system equivalent to equations 5-8:
−k1 0 0
k1 G 0
0 k2 0
−εy2y3
−εy2y3
 (9)
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:
ẋ = Ax+ g(x, y) (10)
ẏ = By + j(x, y) (11)
with (x, y) ∈ Rn × Rm, the n × n matrix A having eigenvalues with zero real
part and the m×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 (9) above is not in normal form since the y′3 equation corresponds
to the zero eigenvalue piece and the nonlinear term of y′3 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
 = T
0 0 G+ k1
0 1 −k1
 . (12)
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:
u′ = 0 · u+
f(u, v, w) (13)
0 −k1
f(u, v, w)
where f(u, v, w) = − ε
(v−k1w)(k2v+G(u+k2w)). 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 6= 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
u′ = 0. (15)
CENTER MANIFOLDS AND LIE SYMMETRIES 4
For the case G = 0, we have a slightly simpler system of equations:
−k1 0 0
k1 0 0
0 k2 0
−εy2y3
−εy2y3
 (16)
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, k2
). Then to transform our
system we again let T be the matrix consisting of these eigenvectors and let
 = T
0 0 1
0 1 −1
1 0 k2
 . (17)
As above, this allows us to write our system in normal form:
g(u, v, w)
g(u, v, w)
w′ = −k1w + 0 (19)
where g(u, v, w) = −ε(v−w)(u+ k2
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 manifoldW 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
ẋ = Ax+ g(x, y) (20)
ẏ = By + j(x, y) (21)
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:
(Mh)(x) = Dh(x)[Ax + g(x, h(x))] −Bh(x)− j(x, h(x)) = 0.
The sign of G does not change the outcome of this calculation, thus we treat
the case G 6= 0 at once. We have h : V → R2, V ⊂ R a neighborhood of the
origin. Thus let h(x) = (h1(x), h2(x)) = (ax
2+ bx3+O(x4), cx2+dx3+O(x4)).
Then f(x, h1(x), h2(x)) = ε(−a+ k1c)x
3 +O(x4) resulting in
(Mh)(x) =
−Gax2 + (−Gb− εa+ k1εc)x
3 +O(x4)
2 + k1dx
3 +O(x4)
CENTER MANIFOLDS AND LIE SYMMETRIES 5
Solving for (Mh)(x) = 0, h1(x) = h2(x) = O(x
4). Thus up to third order, we
have h1(x) = h2(x) = 0, so a center manifold is simply the u−axis.
Next we consider the case G = 0. Here h : V → R, V ⊂ R2, a neighborhood
of the origin. We let h(x) = h(x1, x2) = ax
2+cx1x2+dx
2. Then we calculate
(Mh)(x) =
hx1(x1, x2), hx2(x1, x2)
k2x2 + g(x1, x2, h(x1, x2))
g(x1, x2, h(x1, x2))
+ k1h(x1, x2)
= (2ak2 + ck1)x1x2 + (ak1)x
1 + (ck2 + bk1)x
2 + (dk1)x
1 + (jk2 + ek1)x
+(3dk2 − 2aε− cε+ fk1)x
1x2 + (2fk2 − cε− 2bε+ jk1)x1x
resulting in h(x1, x2) = O(x
4), thus h(x1, x2) = 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
k = ωk(t, y1, y2, . . . , yn), k = 1, . . . , n (22)
the Lie symmetries that transform the variables t, y1, . . . , yn have infinitesimal
generators of the form
X = ξ∂t + η1∂y1 + η2∂y2 + · · ·+ ηn∂yn (23)
where ξ = ξ(t, y1, y2, . . . , yn) and ηk = ηk(t, y1, y2, . . . , yn) for all k. The in-
finitesimal generator must satisfy the Linearized Symmetry Condition:
(1)(y′k − ωk) = 0, k = 1, . . . , n (24)
when (22) holds. In this case the prolongation of X is as follows:
X(1) = X + η
1 ∂y′1 + η
2 ∂y′2 + · · ·+ η
n ∂y′n (25)
where η
is defined as η
= Dtηk − y
kDtξ. The total derivative Dt in this
case is Dt = ∂t + y
1∂y1 + · · ·+ y
n∂yn . Thus we have the following:
= ∂tηk+y
1∂y1ηk+y
2∂y2ηk+· · ·+y
n∂ynηk−y
k(∂tξ+y
1∂y1ξ+y
2∂y2ξ+· · ·+y
n∂ynξ).
A system of first order ODEs has an infinite number of symmetries. We
find symmetries by solving for the functions ξ, ηk that satisfy the Linearized
Symmetry Condition (24). 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
CENTER MANIFOLDS AND LIE SYMMETRIES 6
for our system. We are limited in the symmetries we can calculate by our
computing power. In the case G 6= 0 we calculate symmetries up to third order
in our original coordinates y1, y2, and y3 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 6= 0 case since we want to be able to
compare cases.
4.1 The case G 6= 0
The infinitesimals of the Lie symmetries (up to order 3) are listed in an array
with {ξ, η1, η2, η3}, representing the infinitesimal generator X = ξ∂t + η1∂y1 +
η2∂y2 + η3∂y3 .
X1 = {1, 0, 0, 0}
X2 = {y2,−k1y1y2, k1y1y2 +Gy
2 − εy
2y3, k2y
2 − εy
X3 = {y3,−k1y1y3, k1y1y3 +Gy2y3 − εy2y
3 , k2y2y3 − εy2y
X4 = {0,−y1, y1 +
y2y3,
y2y3}
X5 = {
t,−ty1, ty1 +
ty2 −
ty2y3,
ty2 −
ty2y3}
X6 = {
y1y2 + y1y2y3,
y1y2 + y1y2y3}
Then we transform the infinitesimal generators of the Lie symmetries found in
the yi coordinates as follows. If X is an infinitesimal generator in yi, then X̃ =
(Xt)∂t+(Xu)∂u+(Xv)∂v+(Xw)∂w is the corresponding infinitesimal generator
for a Lie symmetry in the u, v, w coordinates The transformed symmetries in
the form X̃ = {ξ̃, η̃1, η̃2, η̃3} where X̃ = ξ̃∂t + η̃1∂u + η̃2∂v + η̃3∂w:
X̃1 = {1, 0, 0, 0}
X̃2 = {j(u, v, w),
(G− k2)j(u, v, w)f(u, v, w), j(u, v, w)(Gv + f(u, v, w)),−k1wj(u, v, w)}
X̃3 = {l(u, v, w),
(G− k2)l(u, v, w)f(u, v, w), l(u, v, w)(Gv + f(u, v, w)),−k1wl(u, v, w)}
X̃4 = {0,
(G− k2)f(u, v, w),
(Gv + f(u, v, w)),−w}
X̃5 = {
(G− k2)f(u, v, w),
(Gv + f(u, v, w)),−tw}
X̃6 = {m(u, v, w),
(G− k2)m(u, v, w)f(u, v, w),m(u, v, w)(Gv + f(u, v, w)),−k1wm(u, v, w)
where f(u, v, w) is as above, j(u, v, w) = v − k1w, l(u, v, w) = u +
v + k2w
and m(u, v, w) = − 1
(G+ k1)w.
CENTER MANIFOLDS AND LIE SYMMETRIES 7
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 {y1, y2, y3} and list them as X = {ξ, η1, η2, η3},
representing the infinitesimal generator X = ξ∂t + η1∂y1 + η2∂y2 + η3∂y3 .
X1 = {1, 0, 0, 0}
X2 = {y2,−k1y1y2, k1y1y2 − εy
2y3, k2y
2 − εy
X3 = {y3,−k1y1y3, k1y1y3 − εy2y
3 , k2y2y3 − εy2y
X4 = {0,
y1 + y2y3,
y2 + y2y3}
X5 = {
ty1 + ty2y3,
ty2 + ty2y3}
X6 = {
1 + y1y2y3,
y1y2 + y1y2y3}
Then we transform these to the u, v, w coordinate system as above with X̃ =
{ξ̃, η̃1, η̃2, η̃3} where X̃ = ξ̃∂t + η̃1∂u + η̃2∂v + η̃3∂w:
X̃1 = {1, 0, 0, 0}
X̃2 = {v − w,
(v − w)n(u, v, w),−ε(v − w)p(u, v, w),−k1(v − w)w}
X̃3 = {u+
w)n(u, v, w),−ε(u +
w)p(u, v, w),−k1(u+
X̃4 = {0,
n(u, v, w), p(u, v, w),
X̃5 = {
tn(u, v, w), tp(u, v, w),
X̃6 = {
wn(u, v, w), wp(u, v, w),
where n(u, v, w) = εk2w(−v + w) + k1(k2v + εu(−v + w)) and p(u, v, w) =
(v − w)(u + k2
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 ω(u) will be invariant under a given
Lie symmetry, in their notation X = φ∂u + ψ∂v.
CENTER MANIFOLDS AND LIE SYMMETRIES 8
Lemma 5.1 (Lemma 4 of [2] Chapter 7). A center manifold w(u) is invariant
under a Lie symmetry X = φ∂u + ψ∂v if and only if
ψ(u, ω(u)) = (∂u(ω(u))) · φ(u, ω(u)).
For the case G 6= 0, ω(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 X2 with φ(u, v, w) =
(G − k2)j(u, v, w)f(u, v, w) and
ψ(u, v, w) = {j(u, v, w)(Gv + f(u, v, w)),−k1wj(u, v, w)}. Since j(u, 0, 0) ≡ 0,
ψ(u, ω(u)) = ψ(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, ω(u) = 0, and again the right side of our
equation is zero. Thus we must have ψ(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
u′ = k2v − εuv (27)
v′ = −εuv. (28)
The nontrivial symmetries inherited by this system are
X̂2 = {v,
vn(u, v, 0),−εvp(u, v, 0), 0}
X̂3 = {u,
un(u, v, 0),−εup(u, v, w), 0}
X̂4 = {0,
n(u, v, 0), p(u, v, 0), 0}
X̂5 = {
tn(u, v, 0), tp(u, v, 0), 0}
If we transform back to our original variables, we see that on the center
manifold u = y3 and v = y2, resulting in the system:
y′2 = −εy2y3 (29)
y′3 = k2y2 − εy2y3 (30)
CENTER MANIFOLDS AND LIE SYMMETRIES 9
and the symmetries:
X̂2 = y2∂t +
2 − εy
∂y2 − εy
2y3∂y3 (31)
X̂3 = y3∂t +
k2y2y3 − εy2y
∂y2 − εy2y
3∂y3 (32)
X̂4 =
y2 + y2y3
∂y2 + y2y3∂y3 (33)
X̂5 =
t∂t +
ty2 + ty2y3
∂y2 + ty2y3∂y3 (34)
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 X̂4 with η2(t, y2, y3) =
y2 + y2y3 and η3 =
y2y3, this means that, letting γ be the parameter of the one-parameter Lie
group, we need to solve the following for ŷ2 and ŷ3, giving us the map (ŷ2, ŷ3)
as our symmetry:
ŷ2 + ŷ2ŷ3 (35)
= ŷ2ŷ3 (36)
This is equivalent to the system above. Attempting to solve this system we find
it equivalent to solving the following:
ŷ2 = e
+ŷ3)dγ (37)
ŷ3 = e
ŷ2dγ (38)
with the initial conditions ŷ2(γ, y2, y3)|γ=0 = y2 and ŷ3(γ, y2, y3)|γ=0 = y3.
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
ki. In particular, with k3 = 0 and G > 0, they found unrestrained growth of
both M∗ and A. For the values k3 = 0 and a negative G, A increases toward an
upper limit and M∗ increases slightly but then begins to decrease toward zero.
CENTER MANIFOLDS AND LIE SYMMETRIES 10
For k3 > 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 k3 and positive
values of G.
In our consideration of the system, we also found thatM∗(= y3) and A(= y2)
were the two variables that determined the behavior of the system. In the G 6= 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 27 and 28. 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.
k2 = k4 = 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.
References
[1] B. J. Cantwell, Introduction to Symmetry Analysis, Cambridge Univer-
sity Press, Cambridge, United Kingdom, 2002.
[2] G. Cicogna and G. Gaeta, Symmetry and Perturbation Theory in Non-
linear Dynamics, Springer-Verlag, 1999.
[3] E. Ross, I. Taub, C. Doona, F. Feeherry, K. Kustin, The math-
ematical properties of the quasi-chemical model for microorganism growth
– death kinetics in food, International Journal of Food Microbiology, 99
(2005), pp. 157–171.
[4] I. A. Taub, F. E. Feeherry, E. W. Ross, K. Kustin, and
C. J.Doona, A Quasi-Chemical Kinetics Model for the Growth and Death
of Staphylococcus aureus in Intermediate Moisture Bread, Journal of Food
Science, 68, No. 8 (2003), pp. 2530–2537.
	Introduction
	The model
	A simplification
	Normalizing the system
	Center manifold calculations
	Lie Symmetry
	The case G=0
	The case G=0
	The connections between the center manifold and the Lie symmetry
	Comparison to previous results
