One of the difficulties of modeling chemical reactions with differential equation would be that setting the governing equation is not always simple and intuitive. I would suggest you to do some further readings if you are interested in modeling chemical reactions more seriously.

• What governing equation does a reaction follow ?
• How does a reaction sequence become a set of equations ?
• How do the concentrations change over time ?
• Reference

What governing equation does a reaction follow ?

Following table would help you with figuring out a governing equation that will be used for the example in this page. This table is from the reference linked here. I would strongly suggest you to read this reference because it would give you not only for setting up the governing equation but also for analysing a sequence of related chemical reactions at a system level.

The table is the law of mass action. The rate of a reaction is proportional to the product of the concentrations of its reactants, and the rate constant k is the proportionality factor. The first row is a unary reaction, where one molecule of A changes by itself, so [A]' = -k[A]. The second row is a hetero binary reaction, where A₁ meets A₂. Its rate depends on both concentrations, so [A_i]' = -k[A₁][A₂]. The third row is a homeo binary reaction, where two molecules of A meet. Each reaction event uses up two molecules of A, and this is where the factor 2 in -2k[A]² comes from.

The units column tells the same story from a different side. A unary rate constant has the unit s^-1, because k[A] must be a concentration per second. A binary rate constant has the unit M^-1s^-1, because it multiplies two concentrations. The condition at the bottom of the table says that no reactant appears again among the products. When one does, as in A + X -> 2X below, only the net change of that species counts.

• A reaction rate is a product of concentrations : One reactant gives k[A], and two reactants give k[A₁][A₂].
• Stoichiometry sets the coefficient : A reaction that uses two molecules of a species removes that species at twice the reaction rate.
• Count the net change for each species : A species that appears on both sides of a reaction changes only by the difference of its two coefficients.

NOTE : Before look into the derivation of the equation, it would be good to have some intuitive understandings on the solution of the differential equation for this model. Check out this page and see how the solution graph changes as the parameters in the equation changes.  

How does a reaction sequence become a set of equations ?

Now let's look into an example. I got this example reaction sequence from the reference linked here and explained it in my way in this note.

Assume that you have a tank which contains the chemicals named A, X, Y, B and all of these are undergoing chemical reactions as shown below.

Your job is to build a set of differential equations predicting the concentration of each chemicals along with time.

First, let's build a differential equation for the chemical A. To do this, first identify all the chemical reactions which either consumes or produce the chemical (i.e, identify all the chemical reactions in which the chemical A is involved). And then build a differential equation according to the governing equation as shown below.

Next, let's build a differential equation for the chemical X. To do this, first identify all the chemical reactions which either consumes or produce the chemical (i.e, identify all the chemical reactions in which the chemical X is involved). And then build a differential equation according to the governing equation as shown below.

Next, let's build a differential equation for the chemical Y. To do this, first identify all the chemical reactions which either consumes or produce the chemical (i.e, identify all the chemical reactions in which the chemical Y is involved). And then build a differential equation according to the governing equation as shown below.

Next, let's build a differential equation for the chemical B. To do this, first identify all the chemical reactions which either consumes or produce the chemical (i.e, identify all the chemical reactions in which the chemical B is involved). And then build a differential equation according to the governing equation as shown below.

Now we got all the differential equations representing the concentration of each individual chemicals. The last step is to put all the equations together as a simulteneous equation as shown below.

Let's check the system before solving it. If you add the four equations, every term cancels: -k₁AX + (k₁AX - k₂XY) + (k₂XY - k₃Y) + k₃Y = 0. So A + X + Y + B stays constant over time. This is expected, because each reaction turns one molecule into one other molecule. For example, A + X -> 2X turns one A into one X. A conservation check like this finds a wrong sign in any single equation.

X appears on both sides of A + X -> 2X. It enters with the coefficient 1 and leaves with the coefficient 2. So the net gain is one X per reaction, and the term k₁AX in dX/dt carries no factor 2.

The step for the chemical B needs one correction. In the governing law diagram for [B], the rate of the reaction producing [B] is drawn as 0, and the rate consuming [B] is drawn as k₃Y. That gives dB/dt = 0 - k₃Y. But the reaction Y -> B produces B, so the two labels are swapped. The rate producing [B] is k₃Y, and the rate consuming [B] is 0. Its last line, dB/dt = k₃Y, is correct, and it is the one used in the system above.

• One equation per chemical : Each chemical collects a plus term for every reaction that produces it and a minus term for every reaction that consumes it.
• The same reaction term appears in several equations : k₂XY is a loss for X and a gain for Y. A flow between two tanks works the same way.
• The system is nonlinear : The products AX and XY mean that the system cannot be written with a constant matrix, so it is usually solved numerically.

NOTE : Even though it is not exactly same equation, I put the graphical solution of similar differential equation on this note in www.slide4math.com.

How do the concentrations change over time ?

The system above has no closed form solution in general, so let's solve it numerically and read the result. I use k₁ = 1, k₂ = 1 and k₃ = 0.5, with the initial concentrations A = 1, X = 0.01, Y = 0.01 and B = 0. These values are my own choice for illustration. They are not taken from the reference.

The table below lists the four concentrations at selected times, computed with a numerical ODE solver. The last column is the total A + X + Y + B. It stays at 1.02, as the conservation check predicts.

The result runs in three stages. In the first stage, up to about t = 8, X grows at the cost of A. X starts small, but the reaction A + X -> 2X makes more X from each X, so the growth becomes fast once X is large. X reaches its peak of about 0.98 at t = 9.2, when A is almost used up. In the second stage, Y grows from X through X + Y -> 2Y, and it peaks at about 0.159 near t = 17.9. In the third stage, Y turns into B, and B approaches its final value.

Notice that X does not go to zero. The equation for Y can be written as dY/dt = (k₂X - k₃)Y. So Y grows only while X is larger than k₃/k₂ = 0.5. After X falls below 0.5, Y decays, and the conversion of X slows down and stops. About 0.20 of X is left at the end, and B ends near 0.82.

• An autocatalytic step starts slowly and then grows fast : A + X -> 2X uses X to make more X, so its rate grows with X itself.
• A threshold decides whether Y can grow : Y increases only while k₂X > k₃. This threshold is the reason why part of X is never converted.
• Use the conserved total as a numerical check : A numerical solution that does not keep A + X + Y + B constant has a wrong equation or too large a time step.

Reference

[1] From Processes to ODEs by Chemistry (Luca Cardelli, Microsoft Research)

[2] Math 231:Introduction to Ordinary Differential Equations (H. Finotti, The University of Tennessee)

[3] 14.3: Concentration and Rates (Differential Rate Laws)

[4] Chemical reaction rates