13 Simultaneous equations
Contents
13 Simultaneous equations#
13.1 Sequential chemical reactions \(\displaystyle A \stackrel{k_1} \longrightarrow B \stackrel{k_2}\longrightarrow C\)#
Complex chemical reactions can often be represented as a set of simultaneous reactions. The sequential scheme \(\displaystyle A \stackrel{k_1} \longrightarrow B \stackrel{k_2}\longrightarrow C\) has already been solved with the integrating factor method in Section 2, and as an eigenvalue - eigenvector equation in Chapter 7.12.3. Here it is converted into a second-order equation and solved using the \(D\) operator method. The rate equations are,
Differentiating the second equation and substituting \(dA/dt\) produces
Isolating \(A\), substituting and rearranging gives
The solution can be written down after solving the characteristic equation
which has roots \(m = -k_1, -k_2\). The homogeneous equation is
and the constants \(c_1\) and \(c_2\) are, as usual, determined by the initial conditions. If for instance, \(A_0\) is the initial concentration of A, and B and C are initially zero, then since C does not decay, \(C = A_0 - A - B\). Because \(B\) is initially zero, then
and the final step is to find \(c_1\) when
This is done by differentiating B and equating this with the rate equation
after substituting for A and B. The resulting equation produces \(c_1 = k_1A_0/(k_2 - k_1)\) making
which is the same result as that from the integrating factor method.
13.2 Sequential reactions with equal rate constants#
There are important chemical examples of the case when the rate constants are all equal at \(k\) in the scheme \(A \to B \to C\), notably the mechanical unfolding of proteins and DNA using an atomic force microscope. The fact that all the rate constants are equal does not mean that the concentration of B is zero at all times, even though it is formed at the same rate as it decays. Taking the last result, it is not possible to make \(k_1 = k_2\) because the concentration of B becomes undefined as 0/0, however, l’Hopital’s rule could, in this case, be used to find the limit \(k_2 - k_1 \to 0\), (see Chapter 3.8). Instead, we start with the combined rate equations,
The roots of the auxiliary equation are \(m = -k, -k\), and, as they are the same, the method of Section 4.3(iii) has to be used. The first part of the solution is \(\displaystyle B = c_1e^{- kt} + B_1\) where \(B_1\) is a new function of \(t\), which has yet to be found. This solution is normally guessed (from experience), and then tried to see if it fits the equation. As a test, let \(\displaystyle B_1 = c_2te^{- kt}\) making the solution
If B is differentiated and put back into (37), the result is zero showing that this is indeed a solution. To evaluate the constants, this result should be differentiated and \(A\) and \(B\) substituted into the rate equation describing \(B\), as was done in the previous calculation. The initial conditions are normally \(A=A_0\) and \(B = 0\) at \(t = 0\). As \(B = 0\), then \(c_1 = 0\) and, after substituting into the rate equation \(c_2 = kA_0\), this makes
Plotting the A and B concentrations shows that as A falls exponentially, B rises and falls as intuition dictates.
If there are a series of reactions, \(A \to B \to C \to D \to \) etc. and only A is present initially and \(k\) is the same rate constant for each step, then the concentration of the \(n\)th species \(C_n\) is found to be
which is a Poisson distribution.
fig 15a. Concentration of species \(A \to B \to C \to D \to \) etc. calculated using \(\displaystyle C_n =A_0\frac{k^nt^n}{n!} e^{-kt}\) calculated with \(k = 1\) and \(n=0\cdots 5\). At \(t = 0,\, A_0 = 1\) and all other species are zero.
(i) Atomic force microscope unfolding a concatenated protein#
Although rate constants for the reaction of different molecules can be the same, they are generally different. One situation where rate constants are all the same is the unfolding of concatenated proteins by mechanical force.
A concatamer is a macroscopic molecule made by attaching several identical folded proteins in tandem, just like beads on a necklace. The links between proteins consist of only a few peptides compared the the protein which may contain a 100 amino acids. When one end of the concatamer is anchored to a substrate and the other to the tip of an atomic force microscope (AFM), the concatamer can be stretched by the AFM and each protein will unfold but in a random order with respect to its position in the concatamer (Reif et al. 1977; Brockwell et al. 2003). See Fig. 3.21 for a sketch of the experiment. The unfolding is registered by measuring the force and, to a good approximation, the unfolding rate constants for each protein are all the same. However, in this experiment, unlike normal chemical kinetics, the probability or chance \(S\) that a protein is still folded at any given force is measured. If there are four proteins the reaction sequence is \(\displaystyle S_4 \to S_3 \to S_2 \to S_1\) where \(S_4\) is the chance that all proteins are folded, \(S_3\) that three remain folded etc. at a given force.
The rate constants \(k_0\) of thermally induced chemical reactions are constant at fixed temperature and are the same for each protein because they are identical. In the mechanical pulling experiment, the rate constant depends on the force \(f\), because the barrier to unfolding is lowered by the applied force (Bell 1978; Evans & Richie 1997). The rate constant also depends on how fast the force is applied, and the protein becomes stiffer, i.e. more force is needed to unfold it the faster the force is applied. This occurs because the barrier can be lowered more quickly than the average time between thermally induced protein fluctuations that lead to barrier crossing.
The unfolding rate constant is given by
where \(f_0 =kT/x_u\) and \(L\) is the load rate in pN s\(^{-1}\). The constant \(k_0\) is the (thermal) unfolding rate constant at zero force and \(x_u\) is the distance from the minimum of the potential well to the top of the barrier separating the folded and unfolded protein. It is a measure of the distance the protein conformation has to change to reach the unfolding transition state and is typically \(0.25\) nm. The unfolding forces are \(100 \to 300\) piconewton and \(f_0 \approx 15\) pN. Since all proteins are identical, and we assume that each protein experiences the same force irrespective of how many are folded or unfolded, the unfolding scheme for \(N\) proteins can be written as
where the expressions are derivatives with force not time. This set of equations also applies if a set of parallel bonds is unzipped (or unpeeled); they could be the hydrogen bonds holding a protein \(\beta\)-sheet or double stranded DNA together when an opposite force is applied to the top of each strand Fig. 15b.
The set of rather formidable looking simultaneous equations is not linear in \(f\); however, the substitution
greatly simplifies them. The derivative is
and substituting using \(\displaystyle \frac{dS}{ df}\frac{ df}{ du}\) produces,
This set of equations can be solved by an extension of the method of equation (37) and then \(u\) substituted for \(f\). If the initial probability of unfolding is \(S_0\) and all other \(S\) are zero when \(f\) = 0, then
and by continuing the calculation it is found that
and \(S_{unfold}\) is the integral of this result from \(0 \to u\). This equation has the same form as used in figure 15a, but where force replaces time and chance (probability) of unfolding replaces concentration.
Fig. 15b Highly schematic sketch of several parallel hydrogen bonds holding two strands of protein or two lengths of DNA together and being unzipped by an applied force.
13.3 Bloch equations and NMR#
In an NMR experiment, the sample is placed in a large permanent magnetic field aligned along the \(z\)-axis. The magnetic field splits the energy of the relevant nuclear spin states and if circularly polarized RF radiation of the correct frequency is applied this will cause transitions between the nuclear spin energy levels. The time and spatial evolution of the magnetization (the vectorial sum of the spin magnetic moments/unit volume) is described by the Bloch equations, which are basic to NMR (Flygare 1978; Günther 1992; Levitt 2001). The magnetization \(M\) has components in the \(x, y\), and \(z\) directions and decays with a lifetime of \(T_1\), due to longitudinal or population relaxation and by \(T_2\) transverse, or spin-spin, relaxation due to loss of spin coherence. See figure 13 for a sketch of the geometry of an experiment.
The equilibrium magnetization is \(M_0\) and the magnetic field can have components in the \(x,y\), and \(z\) directions, \(B_{x,y,z}\). The Bloch equations describing the experiment are,
Consider now the case where the external field only has a \(z\) component then the equilibrium magnetization vector \(M_0\) points along the z-axis. In a normal NMR experiment, the magnetization is moved from the equilibrium position by the magnetic component of an RF pulse applied in the x - y plane. After the RF pulse has finished (defined as \(t = 0\)) the spins precess only in the external \(B_z\) field starting with the magnetization where it ended up after the RF pulse ended. This will have components along the three axes of \(M_x(0),\; M_y(0)\), and \(M_z(0)\). Because only \(B_z\) is not zero, the Bloch equations are now simplified to
The \(z\) component can be integrated directly with initial condition \(M_z = M_z(0)\) at \(t = 0\) to give
This equation shows that the \(z\) component decays only with a \(T_1\) lifetime, which is the decay of the spin population back to equilibrium. The \(M_x\) and \(M_y\) equations can be integrated as a pair as in section (i), but in this case it is simpler if they are combined as \( M_{+} = M_x + iM_y \) where \(i=\sqrt{-1}\) and afterwards separated into real and imaginary parts. They become
which is a first-order equation and integrates to
and where Euler’s relationship was used to convert the exponentials to sine and cosines. If this equation is expanded and split into real and imaginary parts, \(M_x\) and \(M_y\) are obtained as
where the substitution \(\omega = \gamma B_z\) is also made. This is the NMR transition and Larmor frequency, which is the frequency of precession about the z-axis.
The motion of the magnetization is the vector of the three components. This spirals around the z-axis at the Larmor frequency \(\omega\) until it reaches the equilibrium magnetization pointing along the z-axis, see Fig. 16. The NMR signal is the magnetization’s component as it crosses the x- or y-axis, whichever contains the detecting coil, and this produces the free induction decay or FID, which decays away with lifetime \(T_2\).
To analyse the complete NMR experiment, the effect of the weak field from the coil in the \(x - y\) plane has to be included. In this case, \(B_x\) and \(B_y\) vary as a cosine and sine respectively, and the Bloch equations are far harder to solve. It usual to make the transformation into the rotating frame which means that the axes are changed so that they rotate at the frequency of the R.F. radiation. Flygare (1978), Allen & Eberly (1987), and Günther (1992) discuss this in detail.
The magnetization in a real NMR experiment, while having the form shown in the figure (Fig 16), has many more rotations than could possibly be shown. The typical resonance frequency for a proton is \(400\) MHz and \(T_1\) lifetime \(10\) s. To make the diagrams a frequency of \(0.1\) was used and \(T_1 = 200\) and \(T_2 = 400\) was used in the left figure. The initial position of the vector \(M{(x,y}(0)\) must be defined also. \(M_0\) was calculated as the length of the initial vector. In the figure the initial magnetization vector was at \((-1, -1 , 0)\). The shape of the FID can be seen by plotting \(M_x\) or \(M_y\).
Fig. 16 Magnetization vs. time when \(T_1 = T_2/2 = 200\) (left) and when \(T_1 = 2T_2 = 400\). The magnetization starts at at \(x,\,y,\,z = 1, 1, 0\) which is at the end of the red line, (green dot). This line represents the initial magnetisation. The magnetisation vector ends on the z-axis at point \(M_0 = 1\), red dot.
13.4 A scheme with an equilibrium and reaction . \(A\to C, A\rightleftharpoons B \to C\)#
Reactions where two species could form an equilibrium and at the same time one or the other or both react to form another species are not uncommon. The scheme is related to the ‘circular’ scheme given in chapter 7, 13.2(v).
Our example is of delayed fluorescence of which there are two types. First the inter-molecular process of triplet-triplet annihilation which can produce a spin-statistically limited excited singlet state after two excited triplet states collide (Chapter 10, Q21) and the intra-molecular process of thermally repopulating the singlet state from the triplet and is called thermally activated delayed fluorescence or TADF. This latter process must have a triplet state only slighter lower in energy than the singlet if reverse crossing to the singlet, which is the process that produces delayed fluorescence, is to occur, see scheme fig 16a. The rate constant repopulating the singlet from the triplet is temperature dependent as \(\displaystyle k_r=k_r'e^{-\Delta E/k_BT}\), there \(\Delta E\) is the energy gap, \(k_B\) the Boltzmann constant and \(T\) temperature in degrees Kelvin, see figure 16a. This rate constant and therefore the delayed fluorescence is highly temperature sensitive. The TAPD process is nowadays technologically very important as a process by which to improve the emission efficiency of molecules used as organic light emitting diodes OLEDs although the effect was observed decades earlier in the dye eosin by Parker and Hatchard (Trans. Faraday Soc. 1894, v 57, 1961)) it has until now remained somewhat of a curiosity.
The scheme in figure 16a shows the processes involved and these produce two coupled first-order differential equations. We know that fluorescence with rate constant \(k_f\) is a spin allowed process but phosphorescence \(k_p\) is not, hence \(k_f \gg k_p\). As the reverse intersystem crossing \(k_r\) is thermally activated then \(k_i \gt k_r\). The fluorescence rate constant \(k_f\) and intersystemn crossing \(k_i\) are usually of similar magnitude. The prompt fluorescence decay lifetime is \(\displaystyle \tau_f =\frac{1}{k_f+k_i}\) and the prompt fluorescence yield \(\displaystyle \phi_f=\frac{k_f}{k_f+k_i} =k_f\tau_f\).
figure 16a. Scheme showing the rate constants in thermally activated delayed fluorescence. The term \(\Delta E/k_BT\) needs to be as small as possible to make the rate constant \(k_r\) large enough to significantly repopulate the singlet.
The equations describing the the singlet \(S\) and triplet \(T\) populations are
which we simplify by letting \(k_0 = k_f+k_i\) and \(k_1=k_r+k_p\),
The method to solve these is quite general using the ‘D’ operator method. There are two simultaneous equations which allow us to isolate S or T resulting in a second order equation and then we use the ‘D’ method to find the solution. Rewriting the equations using the ‘D’ operator gives,
It is necessary now to obtain an equation in either S or T. Multiplying the second of equations 37a by \((D+k_0)/k_i\) gives
and subtracting the first equation produces,
which is now a second order differential equation that can be solved for \(T\), and \(S\) can be found by substitution using this solution. To solve this quadratic and only for clarity we let \(a=1,b=k_0+k_i\) and \(c=k_ik_0-k_ik_r\) then
which gives the solution
where \(c_1,c_2\) are constants to be determined from the initial conditions. These are that at \(t=0, S=S_0, T=0\) thus \(c_1=-c_2\) and therefore substituting gives
If we differentiate this solution then we can use \(\displaystyle dT/dt=-k_iS+k_1T\) to find \(S\). Doing this gives
and from which we can find \(c_1\) because \(S=S_0\) when \(t=0\) making \(\displaystyle c_1=k_iS_0/q\) and finally
After some very fiddly algebra,
where \(q\) is given above.
We know from our statement of the kinetics that the triplet is not present initially and must decay to zero at long times, and we can see from the form of these equations that the triplet population is indeed zero at \(t=0\) and zero at \(t=\infty\) so must reach a maximum at some finite time and so is described by the difference of two exponentials. The singlet decays from its initial value \(S_0\) with an initial decay due to ‘prompt’ fluorescence and a longer decay due to re-population of \(S_1\) from the triplet, causing the delayed fluorescence, so the singlet is expected to decay as the sum of two exponentials, as is confirmed by its equation. The two excited state lifetimes are
and as \(q\) is positive \(\tau_1\) this must be the shorter lifetime. Using the values to generate figure 16b, \(k_f=0.1, k_i=0.1\), \( k_r = 0.1e^{-4} = 1.8\cdot10^{-3},\) \( k_p = 5\cdot10^{-4}\) ns\(^{-1}\) then \(\tau_1 = 4.98\) ns, which is almost the same as \(1/k_0\), and \(\tau_2 = 709\) ns. The large difference in these values is the result of the large energy gap from the triplet to the ground state and spin-forbidden nature of the phosphorescence.
Figure 16b. Prompt fluorescence (solid grey line), total fluorescence (red) and phosphorescence (dashed) vs. time on a log-log scale. The delayed fluorescence is the difference between the prompt and total fluorescence. The rate constants used were (in units of ns\(^{-1}\)), \(k_f=0.1,k_i=0.1,k_r=0.1e^{-4}=1.8\cdot10^{-3}, k_p=5\cdot10^{-4}\)
The direct and delayed fluorescence, and the phosphorescence yields can be calculated by integrating the S and T decay profiles either analytically or numerically. Doing this for the rate constants used in the figure gave a direct fluorescence yield of \(0.5\) and an additional delayed yield of \(0.32\) which is a considerable increase. This may seem strange looking at the figure because the prompt fluorescence has decayed by so much before the delayed emission seems to start. However, the log time-scale is misleading and because this emission is present for such a long time, compared to the prompt fluorescence, and this makes its yield larger than it appears to be. The delayed fluorescence yield can also be worked out analytically, (C. Baleizao & M. Berberan-Santos, J. Chem. Phys, 2007, v. 126, p. 204510). For the purposes of calculation they supposed that the energy moved sequentially from singlet to triplet as \(S^{(1)}\to T^{(1)} \to S^{(2)}\to T^{(2)} \cdots\) and worked out the chance of returning to the singlet at each step \(1,2,3 \cdots\), and so fluorescing at each. The initial chance is simply the fluorescence yield \(\phi_f\) but now a fraction \(\phi_i\) moves to the triplet and a fraction \(\phi_r\) of this returns and so the total chance so far is \(\phi_f(1+\phi_i\phi_r)\). At the next step only the fraction \((\phi_i\phi_r)^2\) returns and by induction \((\phi_i\phi_r)^3\) and so on for further steps, giving
where \(\displaystyle\phi_f=\frac{k_f}{k_f+k_i}\) is the prompt yield and \(\displaystyle\phi_r=\frac{k_r}{k_p+k_r}\) and \(\displaystyle\phi_i=\frac{k_i}{k_f+k_i}=1-\phi_f\). The delayed fluorescence yield can be worked out from the total fluorescence by subtracting the prompt fluorescence and is
To produce a large total fluorescence yield the general conditions are \(k_f \lt k_i, k_r\lt k_f\) and \( k_p\) and \(k_r\) are small than compared to other rate constants. If \(k_r\) can be made comparable to \(k_i\) or larger then the total yield can be increased but this is not generally possible unless unrealistically high temperatures are used.
13.5 Second-order equations. Coupled springs#
Pairs of second-order equations can be solved using the operator method. In this example, the motion of a pair of masses and springs is calculated. This calculation is an alternative to the matrix method of Chapter 7.12.
Consider two identical springs each supporting a mass as shown in Fig. 17. This could be a simple model for part of a vehicle’s suspension. The displacement from equilibrium of the upper mass \(m_1\) is \(y\), and that of the lower one \(x\). The force constants are \(k_1\) and \(k_2\). If the springs were isolated, they would exert a force equal to \(-k_2x\) or \(-k_1y\) on their respective masses. When connected together, the upper spring now exerts a force equal to \(-k_1(y - x)\) on \(m_1\) as the displacement is changed by the lower spring. The total force on the lower spring is also changed and is \(-k_2x + k_1(y - x)\). Together, these produce the force equations,
Fig. 17 Two coupled springs.
To make the algebra simpler, suppose that the masses are each \(m\) and the spring force constants, \(k\). Differentiating twice gives
which reduces to
This can be changed into the \(D\) operator form
where \(\displaystyle \omega^2 = k/m\), and has the characteristic equation
This has four solutions
These solutions are equivalent to
and produces the homogeneous equation
The constants are determined by the initial conditions, which are the initial position, initial velocity, acceleration, and its derivative. The result is that the motion of the spring is the sum of two simple harmonic motions of different frequencies, one greater and one smaller than \(\omega\), which is the frequency of each isolated single spring. Using the second set of solutions, the frequencies are in the golden ratio and its reciprocal, with respect to \(\omega\),
The motion is shown in the next figure (17a) with initial condition \(x_0 = 1\). The motion is clearly complex but appears to repeat itself after about \(16\) secs but more exactly at about \(26\) secs. The figure on the right shows the phase plane with dots at \(0\) and \(16.6\) s (red dot). The red line shows how the trajectory is similar after approx \(16\) secs but not exactly the same as before. Poincare has postulated that any system will periodically return arbitrarily close to its starting conditions; this is seen to be the case in this simple system. The greater the number of oscillators the longer the time for recurrence to occur.
Fig 17a Left. Motion of displacement \(x\) of two coupled weights (see fig. 17) vs. time. Right. Phase plot for the motion showing that the trajectory almost but not exactly arrives back at the initial displacement after \(\approx 16.6\) s. The red line shows that the subsequent motion differs slightly from that which occurred initially.