Nature exhibits a curious asymmetry between the left and the right, which is generally called ‘chira — Physical Chemistry — Kinetics Chemistry Question

Nature exhibits a curious asymmetry between the left and the right, which is generally called ‘chiral asymmetry’. Indeed, living organisms contain mostly L-amino acids and D-carbohydrates. One of the possible explanations of this phenomenon is based on the idea of autocatalysis. Chiral (asymmetric) autocatalysis is a reaction in which every chiral product serves as the catalyst of its own formation. In such reactions small initial excess of one of the enantiomers can increase exponentially in time.

Consider the kinetic scheme explaining this phenomenon. Two Enantiomers, XL and XD, are reversibly formed from achiral reagents T and S:

S + T <=> XL (1) (with rate constants k₁, k_-1)
S + T <=> XD (2) (with rate constants k₁, k_-1)
S + T + XL <=> 2 XL (3) (with rate constants k₂, k_-2)
S + T + XD <=> 2 XD (4) (with rate constants k₂, k_-2)
XL + XD → P (5) (with rate constant k₃)

Enantiomers react with each other giving the product P. The reactions take place in an open system, where constant concentrations of reagents S and T are maintained.

The system of rate equations can be solved numerically using any of the mathematical packages, for example Mathematica, MathCad, etc. Alternatively, you may use the program KINET posted on the official website www.icho₃₉.chem.msu.ru. Let us assume the following values of rate constants (in arbitrary units): k₁ = 0.5, k–1 = 0.1, k₂ = 0.5, k–2 = 0.2, k₃ = 0.5.

Procedure
For numerical solution of the systems of differential equations mathematical packages use different commands. In Mathematica it is done by the function NDSolve. The arguments are the list of equations, initial conditions and a time interval. For example, the system of equations:
a'(t) = -a(t)p(t)
p'(t) = a(t)p(t) - 2p(t)

with the initial conditions a(0) = 2, p(0) = 0.5 in a time interval from t = 0 to t = 10 is solved numerically by the command:
sol=NDSolve [{a ' ==-a *p , p ' == a *p -2*p , a == 2, p == 0.5}, {a, p}, {t, 0,10}]

The obtained solution is presented on the graph by the command Plot:
Plot [Evaluate [{a , p }/.sol, {t, 0,10}], PlotRange^-> All]