MATRIX EXPONENTIAL SOLUTIONS OF FIRST-ORDER CHEMICAL NETWORKS

The matrix exponential method as implemented in MATLAB is demonstrated as a facile tool for solving the time-dependent concentrations of an arbitrary chemically reactive network modelled as a coupled linear system of first-order differential equations. The method is used to verify a 10 species network incorporating experimentally supplied forward rate constants; and a random 11 species network incorporating both forward and backward rate constants as modelled by semiclassical electron transfer theory. Also demonstrated is the matrix exponential solution in exact arithmetic (via Putzer’s algorithm and verified by Laplace transforms) for a chain of three species coupled reversibly.


Introduction
Huang [3] solved via computer simulation the time-dependent concentration curves for a dechlorination and elimination reaction involving 10 species. The technique involved transformation of a coupled linear system of first-order differential equations into an uncoupled system via diagonalization of the reaction rate constant matrix, a technique admitted to suffer the shortcoming of the necessity of a full set of eigenvectors. In contrast, matrix exponential does not suffer from this restriction [5], [7]. Moreover, it is easy to implement in MATLAB as e At can be calculated by the command expm(At) and postmultiplication by a column vector of initial concentrations produces the concentration vector at time t for the network (assuming the process starts at t = 0). As pointed out by Huang [3], the existing literature (e.g., [1], [4], [6]) provides theoretical treatments of matrix methods for chemical systems but lacks concrete computer implementations and sample concentration-time curves, which can be beneficial to both chemical engineers and chemists who desire simple and efficient programs dedicated to first-order kinetic analysis. This paper employs the matrix exponential technique, which can be compared to the diagonalization method [3]. Before providing a general network description, we begin with kinetic equations for a system of three species and provide its solution in exact arithmetic.

Reversibly coupled kinetic equations
Consider the following reversibly coupled chain of three species If we let x, y, and z denote the time-dependent concentrations of A, B, and C, the differential equations for first-order kinetics arė In matrix format this reads We pursue an exact solution for the situation when only A (no B or C) is present initially (for comparison, see [8]). Note thatẋ +ẏ +ż = 0 and so by linearity Since our initial value problem has a unique solution the constant is the initial amount of A. Using a unit amount and choosing x and z as independent variables we obtain y = 1 − x − z so we only need to solve for x and z. The equivalent linear system for two concentration variables becomes The matrix exponential e At is a matrix function defined on square matrix A multiplied by indeterminate parameter t, which represents time in our application; it is suited to solve both the homogeneous (1) and inhomogeneous (2) linear systems presented in section 2. Specifically, e At is defined by the following power series where I is the identity matrix. Fundamental to vector space theory is that e At solves (1) and (2) as follows. The homogeneous system (1) using 3 × 3 matrix A and its solution are given bẏ is a column vector function of the same dimension as ⃗ x(t) and A now represents the associated 2 × 2 matrix. The unique solution is Note that for system (2) we have chosen t 0 = 0 and ⃗ f (t) is a constant vector (the variable t does not even appear).

Putzer's algorithm
Although it is not evident from the power series definition, any matrix exponential function can be computed in exact arithmetic as a time-weighted linear combination of finitely many matrices [7]. Putzer's algorithm proves this claim. A description of the method follows (proof is straightfoward but omitted): Let λ 1 , . . . , λ n be the eigenvalues, possibly complex or repeated, of n × n matrix A. Then Putzer's algorithm states that is a vector function whose components satisfy the following initial value problem: 5 Exact solution of 2 × 2 inhomogeneous system It will be a bit too tedious to deal with the 2 × 2 matrix of constants from (2) directly so instead we solve the system⃗ where a, b, c, d, F, and G are real-valued constants and ⃗ The first step is to calculate the eigenvalues of 2×2 matrix A, which is done with the well-known formula We have two cases to consider depending on whether the eigenvalues are distinct or repeated. We begin with the former case.
In the case that (a+d) 2 = 4(ad−bc) equation (5) produces two repeated eigenvalues. Removing the subscript and calling this eigenvalue λ, Putzer's algorithm now gives p 1 = e λt p 2 = te λt and by working through similar steps one can verify the following: A and ⃗ f are given by (4) and λ 1 = λ 2 ≡ λ is given by For computation of the exact solution it is perhaps easiest to implement a function (in say MATLAB or EXCEL) that takes parameters a, b, c, d, F, G and time t as inputs and produces the concentration vector as the output by substituting the appropriate combination of rate constants for the inputs according to kinetic system (2).
The Laplace transform method can be used to obtain the exact solution of (4) and verify the results (6), (7), (8), (9) from Putzer's method. We outline how this is done for z(t) when λ 1 ̸ = λ 2 ; the remaining cases for x(t) and λ 1 = λ 2 are similar. The first step involves taking the Laplace transform of the inhomogeneous system (4). Letting x(t) ↔ X(s) and z(t) ↔ Z(s) be Laplace transform pairs we find By Cramer's rule we can solve for X(s) and Z(s). The result for Z(s) is Partial fractions and inversion give Using the identity λ 1 + λ 2 = a + d we can show that z(t) matches (7), so the methods of Laplace and Putzer agree, as required.

Network modelling
For purpose of solving a general first-order network via matrix exponential it is straightforward to handle the homogeneous system⃗ x = A⃗ x directly and compute ⃗ x(t) = e A(t−t0) ⃗ x 0 . The reaction rate constant matrix A can be found as follows. For a network consisting of only species i and j are coupled reversibly i kij −−⇀ ↽−− kji j we havė For a general network species i may be coupled to other species thus A ii =j k ij where the summation is over all species j coupled to species i. The constant k ji stands alone at row i column j so that A ij = k ji . Likewise A jj =i k ji and A ji = k ij . If we model i and j as irreversibly coupled i kij −−→ j the above considerations apply using k ji = 0. Construction of A in MATLAB can be performed by preconstructing a matrix {k ij } of rate constants. When k ij is encountered in a traversal make the assignments where A is initialized to zero.
7 Solution curves for a 10 species network [3] In this network the species are irreversibly coupled with arrows indicating the reaction direction according to the diagram (Figure 1). Figure 1: Diagram of 10 species network [3]. Rate constants in units h -1 . Species 1 is assigned 1M initial concentration, the other species zero.

11 species network modelled with semiclassical electron transfer theory
In this theory we view the network as an electric circuit in which free energies assigned to the species act as electric potentials and distances assigned to coupled pairs (i.e., straight-line distances between connected species) act as resistances. The semiclassical formula [9] for the rate constants of reversible connection i kij −−⇀ ↽−− kji j is given by (swap i and j to obtain formula for where H is constant (s -1 ), r ij = r ji is the Angstrom distance, ∆G • ij = G • j − G • i is the difference in absolute standard potentials (measured in volts), λ is constant (called the reorganization energy, measured in volts), and k B T is the Boltzmann factor. For this study was used H = 1.5e14s -1 , λ = 0.7V , and T = 300K.
For our purpose we read (*) as saying that the passage of current is reduced exponentially as r ij is increased; and that there is Gaussian dependence of rate on potential difference, with the maximum driving force occuring at ∆G • ij = -λ = -0.7V , below which remarkably the rate begins to slow down rather than speed up (this is the Marcus inverted region effect; see [2]). For this study was chosen an 11 species network whose graph is seen here: Nonzero potentials of 0.1V are assigned to sources 1, 2, and 3 at initial concentrations 0.2, 0.4, and 0.6M, respectively (other species assigned zero initial concentration). Lower potentials -0.3, -0.2, and -0.2V are assigned to 9, 10, and 11. All other species are assigned zero potential (labels omitted). All distances are 8Å except for the four terminal branches at 15Å.