I_abcd
Derivation of differential equations describing evolution of spin concentrations
1. Reaction rates and partial conversion rates
4. Expression in terms of spin (monomer) concentrations
clean up workspace
reset()
Write properly balanced reactions equations:
Transition A1:
(1) (2)
Ra<=>Rb
Constants: k_1_A_1 (forward), k_2_A_1 (reverse).
Transition A2:
(3) (4)
Rc<=>Rd
Constants: k_1_A_2 (forward), k_2_A_2 (reverse).
Transition B1:
(1) (3)
Ra<=>Rc
Constants: k_1_B_1 (forward), k_2_B_1 (reverse).
Transition B2:
(2) (4)
Rb<=>Rd
Constants: k_1_B_2 (forward), k_2_B_2 (reverse).
Write reaction rates
Introduction.
We distinguish reaction rates (Rate, elementary reaction acts per unit time) and conversion rates (dc/dt, number of moles of the specific species consumed/produced per unit time). Conversion rates, dc/dt, for species are related to reaction rates, Rate, through molecularity coefficients.
To compute conversion rates, we need to distinguish partial conversion rates from net (overall) conversion rates. The net conversion rate is actual rate of change in measured concentration of the species. Partial conversion rate is a conversion rate of the species observed along the specific branch of the reaction mechanism. Summing partial conversion rates of the species one obtains the net conversion rate for this species.
A1 Isomerization forward transition (1_A_1)
a reaction rate
eq1_1a:= Rate_1_A_1 = k_1_A_1*Ra
a partial conversion rate of Ra : one reaction act uses one molecule of Ra
eq1_1b:= dcRadt_1_A_1 = Rate_1_A_1 * (-1)
The final form
eq1_1c:= eq1_1b | eq1_1a
a partial conversion rate of Rb : one reaction makes one molecule of Rb
eq1_1d:= dcRbdt_1_A_1 = Rate_1_A_1 * (+1)
The final form
eq1_1e:= eq1_1d | eq1_1a
A1 Isomerization reverse transition (2_A_1)
a reaction rate
eq1_2a:= Rate_2_A_1 = k_2_A_1*Rb
a partial conversion rate of Ra : one reaction act makes one molecule of Ra
eq1_2b:= dcRadt_2_A_1 = Rate_2_A_1 * (+1)
The final form
eq1_2c:= eq1_2b | eq1_2a
a partial conversion rate of Rb : one reaction act uses one molecule of Ra
eq1_2d:= dcRbdt_2_A_1 = Rate_2_A_1 * (-1)
The final form
eq1_2e:= eq1_2d | eq1_2a
A2 Isomerization forward transition (1_A_2)
a reaction rate
eq1_3a:= Rate_1_A_2 = k_1_A_2*Rc
a partial conversion rate of Rc : one reaction act uses one molecule of Rc
eq1_3b:= dcRcdt_1_A_2 = Rate_1_A_2 * (-1)
The final form
eq1_3c:= eq1_3b | eq1_3a
a partial conversion rate of Rd : one reaction makes one molecule of Rd
eq1_3d:= dcRddt_1_A_2 = Rate_1_A_2 * (+1)
The final form
eq1_3e:= eq1_3d | eq1_3a
A2 Isomerization reverse transition (2_A_2)
a reaction rate
eq1_4a:= Rate_2_A_2 = k_2_A_2 * Rd
a partial conversion rate of Rc : one reaction act makes one molecule of Rc
eq1_4b:= dcRcdt_2_A_2 = Rate_2_A_2 * (+1)
The final form
eq1_4c:= eq1_4b | eq1_4a
a partial conversion rate of Rd : one reaction uses one molecule of Rd
eq1_4d:= dcRddt_2_A_2 = Rate_2_A_2 * (-1)
The final form
eq1_4e:= eq1_4d | eq1_4a
B1 Isomerization forward transition (1_B_1)
a reaction rate
eq1_5a:= Rate_1_B_1 = k_1_B_1 * Ra
a partial conversion rate of Ra : one reaction act uses one molecule of Ra
eq1_5b:= dcRadt_1_B_1 = Rate_1_B_1 * (-1)
The final form
eq1_5c:= eq1_5b | eq1_5a
a partial conversion rate of Rc : one reaction act makes one molecule of Rc
eq1_5d:= dcRcdt_1_B_1 = Rate_1_B_1 * (1)
The final form
eq1_5e:= eq1_5d | eq1_5a
B1 Isomerization reverse transition (2_B_1)
a reaction rate
eq1_6a:= Rate_2_B_1 = k_2_B_1 * Rc
a partial conversion rate of Ra : one reaction act makes one molecule of Ra
eq1_6b:= dcRadt_2_B_1 = Rate_2_B_1 * (1)
The final form
eq1_6c:= eq1_6b | eq1_6a
a partial conversion rate of Rc : one reaction act uses one molecule of Rc
eq1_6d:= dcRcdt_2_B_1 = Rate_2_B_1 * (-1)
the final form
eq1_6e:= eq1_6d | eq1_6a
B2 Isomerization forward transition (1_B_2)
a reaction rate
eq1_7a:= Rate_1_B_2 = k_1_B_2 * Rb
a partial conversion rate of Rb: one reaction act uses one molecule of Rb
eq1_7b:= dcRbdt_1_B_2 = Rate_1_B_2 * (-1)
the final form
eq1_7c:= eq1_7b | eq1_7a
a partial conversion rate of Rd : one reaction act makes one molecule of Rd
eq1_7d:= dcRddt_1_B_2 = Rate_1_B_2 * (1)
the final form
eq1_7e:= eq1_7d | eq1_7a
B2 Isomerization reverse transition (2_B_2)
a reaction rate
eq1_8a:= Rate_2_B_2 = k_2_B_2 * Rd
a partial conversion rate of Rb : one reaction act makes one molecule of Rb
eq1_8b:= dcRbdt_2_B_2 = Rate_2_B_2 * (1)
the full form
eq1_8c:= eq1_8b | eq1_8a
a partial conversion rate of Rd : one reaction act uses one molecule of Rd
eq1_8d:= dcRddt_2_B_2 = Rate_2_B_2 * (-1)
the full form
eq1_8e:= eq1_8d | eq1_8a
To define evolution of the species we need to compute concentrations as a function of time. To this end, we will write differential equations for conversion rates of all species.The net conversion rate of the species is a sum of partial conversion rates along all branches.
Net conversion rate of Ra
Sum all pertaining partial conversion rates
eq3_1a:= dcRadt_N = dcRadt_1_A_1 + dcRadt_2_A_1 + dcRadt_1_B_1 + dcRadt_2_B_1
substitute partial conversion rates
eq1_1c;
eq1_2c;
eq1_5c;
eq1_6c;
eq3_1b:= eq3_1a | eq1_1c | eq1_2c | eq1_5c | eq1_6c;
Net conversion rate of Rb
Sum all pertaining partial conversion rates
eq3_2a:= dcRbdt_N = dcRbdt_1_A_1 + dcRbdt_2_A_1 + dcRbdt_1_B_2 + dcRbdt_2_B_2
substitutions
eq1_1e; eq1_2e; eq1_7c; eq1_8c;
eq3_2b:= eq3_2a | eq1_1e | eq1_2e | eq1_7c | eq1_8c
Net conversion rate of Rc
Sum all pertaining partial conversion rates
eq3_3a:= dcRcdt_N = dcRcdt_1_B_1 + dcRcdt_2_B_1 + dcRcdt_1_A_2 + dcRcdt_2_A_2
substitutions
eq1_3c; eq1_4c; eq1_5e; eq1_6e;
eq3_3b:= eq3_3a | eq1_3c | eq1_4c | eq1_5e | eq1_6e
Net conversion rate of Rd
Sum all pertaining partial conversion rates
eq3_4a:= dcRddt_N = dcRddt_1_A_2 + dcRddt_2_A_2 + dcRddt_1_B_2 + dcRddt_2_B_2
substitutions
eq1_3e; eq1_4e; eq1_7e; eq1_8e;
eq3_4b:= eq3_4a | eq1_3e | eq1_4e | eq1_7e | eq1_8e
not needed here because we do not have oligomerization reactions: one spin in each reactant is converted to one spin in the product.
Summarize the derivation results
eq3_1b
eq3_2b
eq3_3b
eq3_4b
Assign order to species
eq5_1a:= Ra = C1;
eq5_1b:= Rb = C2;
eq5_1c:= Rc = C3;
eq5_1d:= Rd = C4;
Same order for the net rates
eq5_2a:= dcRadt_N = dC1dt;
eq5_2b:= dcRbdt_N = dC2dt;
eq5_2c:= dcRcdt_N = dC3dt;
eq5_2d:= dcRddt_N = dC4dt;
Restate the equations in terms of numbered species
eq5_3a:= eq3_1b | eq5_1a | eq5_1b | eq5_1c | eq5_1d | eq5_2a | eq5_2b | eq5_2c | eq5_2d
eq5_3b:= eq3_2b | eq5_1a | eq5_1b | eq5_1c | eq5_1d | eq5_2a | eq5_2b | eq5_2c | eq5_2d
eq5_3c:= eq3_3b | eq5_1a | eq5_1b | eq5_1c | eq5_1d | eq5_2a | eq5_2b | eq5_2c | eq5_2d
eq5_3d:= eq3_4b | eq5_1a | eq5_1b | eq5_1c | eq5_1d | eq5_2a | eq5_2b | eq5_2c | eq5_2d
Prepare results for transfer to MATLAB
To avoid typing errors when transfering derived K matrix to MATLAB we type it in here and then directly test against the derivation result from above. After that the K matrix may be transfered to MATLAB by cut-and-paste of the MuPad output.
Enter the K-matrix looking at the above results (collect terms at correspondingly numbered species).
Simple rules that allow catching mistakes in K matrix derivation:
(1) a sum of each column should be zero (so each constant must appear with both positive and negative sign), and
(2) each row has to have complete pairs of constants (i.e., if k12
appears there must be k21 in the same row with an opposite sign and so on).
K:=matrix(4,4,[
[ (-k_1_A_1-k_1_B_1), k_2_A_1, k_2_B_1, 0 ],
[ k_1_A_1, -k_2_A_1-k_1_B_2, 0, k_2_B_2 ],
[ k_1_B_1, 0, -k_1_A_2-k_2_B_1, k_2_A_2 ],
[ 0, k_1_B_2, k_1_A_2, -k_2_A_2-k_2_B_2 ]
])
![matrix([[- k_1_A_1 - k_1_B_1, k_2_A_1, k_2_B_1, 0], [k_1_A_1, - k_1_B_2 - k_2_A_1, 0, k_2_B_2], [k_1_B_1, 0, - k_1_A_2 - k_2_B_1, k_2_A_2], [0, k_1_B_2, k_1_A_2, - k_2_A_2 - k_2_B_2]])](I_abcd_images/math80.png)
Create a column vector containing concentrations of species in numbered notation
P:=matrix(4,1,[C1, C2, C3, C4])
![matrix([[C1], [C2], [C3], [C4]])](I_abcd_images/math81.png)
Check correctness of the entered K matrix by multiplying with P and comparing to the above equations:
Multiply K and P:
dCdt_manual_input:= K*P
![matrix([[C2*k_2_A_1 + C3*k_2_B_1 - C1*(k_1_A_1 + k_1_B_1)], [C1*k_1_A_1 + C4*k_2_B_2 - C2*(k_1_B_2 + k_2_A_1)], [C1*k_1_B_1 + C4*k_2_A_2 - C3*(k_1_A_2 + k_2_B_1)], [C2*k_1_B_2 + C3*k_1_A_2 - C4*(k_2_A_2 + k_2_B_2)]])](I_abcd_images/math82.png)
Collect right-hand-side parts of equations
dCdt_mupad:=matrix(4,1,[ rhs(eq5_3a), rhs(eq5_3b), rhs(eq5_3c), rhs(eq5_3d)])
![matrix([[C2*k_2_A_1 - C1*k_1_B_1 - C1*k_1_A_1 + C3*k_2_B_1], [C1*k_1_A_1 - C2*k_1_B_2 - C2*k_2_A_1 + C4*k_2_B_2], [C1*k_1_B_1 - C3*k_1_A_2 - C3*k_2_B_1 + C4*k_2_A_2], [C2*k_1_B_2 + C3*k_1_A_2 - C4*k_2_A_2 - C4*k_2_B_2]])](I_abcd_images/math83.png)
Compare the derivation result to manual input
dCdt_mupad=dCdt_manual_input:
normal(%);
bool(%)
![matrix([[C2*k_2_A_1 - C1*k_1_B_1 - C1*k_1_A_1 + C3*k_2_B_1], [C1*k_1_A_1 - C2*k_1_B_2 - C2*k_2_A_1 + C4*k_2_B_2], [C1*k_1_B_1 - C3*k_1_A_2 - C3*k_2_B_1 + C4*k_2_A_2], [C2*k_1_B_2 + C3*k_1_A_2 - C4*k_2_A_2 - C4*k_2_B_2]]) = matrix([[C2*k_2_A_1 - C1*k_1_B_1 - C1*k_1_A_1 + C3*k_2_B_1], [C1*k_1_A_1 - C2*k_1_B_2 - C2*k_2_A_1 + C4*k_2_B_2], [C1*k_1_B_1 - C3*k_1_A_2 - C3*k_2_B_1 + C4*k_2_A_2], [C2*k_1_B_2 + C3*k_1_A_2 - C4*k_2_A_2 - C4*k_2_B_2]])](I_abcd_images/math84.png)
=> If TRUE ---the typed K-matrix is correct.
Use this K-matrix (copy-paste output to MATLAB)
K;
![matrix([[- k_1_A_1 - k_1_B_1, k_2_A_1, k_2_B_1, 0], [k_1_A_1, - k_1_B_2 - k_2_A_1, 0, k_2_B_2], [k_1_B_1, 0, - k_1_A_2 - k_2_B_1, k_2_A_2], [0, k_1_B_2, k_1_A_2, - k_2_A_2 - k_2_B_2]])](I_abcd_images/math86.png)
I derived differential equations governing spin populations. The K matrix has been prepared for transferring to MATLAB.