Predict

The predict module implements the DRL class which helps to create a prediction of a chemical system. When the prediction fails it raises a InvalidPredictionError. The DRL class implements methods to:

  1. Predict the concentrations for a DRL experiment using the ODE solver (preferred).

  2. Predict the concentrations for a DRL experiment using the explicit Euler formula (discouraged).

  3. Calculate the change in chemical concentration as a function of the current concentrations.

  4. Calculate the Jacobian, which is required by the ODE solver. This only works for reaction where each reaction step is first order in each chemical.

The details how the rate equations and Jacobian matrix are calculated can be found in the implementation details section.

class predict.DRL(reactions, rate_constants, output_order=None, verbose=False)

Contains all information required to predict a chemical system’s concentrations over time.

Parameters:

  • reactions (list[tuple[str, list[str], list[str]]]) – A list of each reaction step that describes the total system. Each reaction step is a tuple, where the first element is the name of the rate constant. The second element contains a list with the names of each reactant. The third element contains a list with the names of each product.

  • rate_constants (dict[str, float] | pd.Series) – The rate constants and their respective values.

  • output_order (list[str]) – The index of the initial concentrations and prediction. The order of each chemical must be given. If None (default), the order will be alphabetical.

  • verbose (bool) – If true, it will print store information on the reactions in the model. This information is also stored as the attribute ‘reactions_overview’.

Variables:

reactions_overview (pd.DataFrame) – If verbose was True upon initialization, this will yield an easier to read overview of the reactions in the system. It also shows the value of each rate constant, and not only its name.

predict_concentration(t_eval_pre, t_eval_post, initial_concentrations, labeled_concentration, dilution_factor, atol=1e-10, rtol=1e-10)

Predicts the concentrations during a DRL experiment. It utilizes the ODE solver ‘scipy.integrate.solve_ivp’ with the Radau method.

Parameters:

  • t_eval_pre (np.ndarray) – The time steps before the addition of the labeled compound. The first element will be the starting time, and the last the time when it ends. It can be a 2-cell array.

  • t_eval_post (np.ndarray) – The time steps after the addition of the labeled compound, that must be evaluated.

  • initial_concentrations (dict[str, float]) – The initial concentrations of each chemical. Only non-zero concentrations are required.

  • labeled_concentration (dict[str, float]) – The concentration of the labeled chemical. This concentration is not diluted.

  • dilution_factor (float) – The factor (≤ 1) by which the prediction will be ‘diluted’ when the labeled chemical is added.

  • atol (float) – The absolute tolerances for the ODE solver.

  • rtol (float) – The relative tolerances for the ODE solver.

Returns:

The predicted concentrations for each time stamp in the t_eval_post array. The time array itself will be appended to the DataFrame.

Return type:

pd.DataFrame

predict_concentration_Euler(t_eval_pre, t_eval_post, initial_concentrations, labeled_concentration, dilution_factor, steps_per_step=1)

Predicts the concentrations during a DRL experiment.

Warning

This method is less accurate and slower compared to using an ODE solver such as implemented in predict_concentration().

Parameters:

  • t_eval_pre – The time steps that must be evaluated, before the addition of the labeled compound.

  • t_eval_post – The time steps that must be evaluated and returned, after the addition of the labeled compound.

  • initial_concentrations (dict[str, float]) – The initial concentrations of each chemical. Only non-zero concentrations are required.

  • labeled_concentration (dict[str, float]) – The concentration of the labeled chemical. This concentration is not diluted.

  • dilution_factor (float) – The factor (<= 1) by which the prediction will be ‘diluted’ when the labeled chemical is added.

  • steps_per_step – The number of steps that should be modeled for each point that is evaluated according to the t_eval arrays.

Returns:

The prediction of the concentration as a function of time after the addition of the labeled compound.

Return type:

pd.DataFrame

calculate_step(_, y)

Calculates the rate of change in the chemical system.

Parameters:

  • _ (ndarray) – Time is inputted here by scipy.integrate.solve_ivp, but it is not used to calculate the rate of change.

  • y (ndarray) – The current concentrations of each chemical.

Returns:

The change in concentration with respect to time. This has NOT been multiplied with the change in time yet!

Return type:

np.ndarray

calculate_jac(_, y)

Calculates the Jacobian for the chemical system. This function is required by the stiff ODE solvers, such as Radau, in scipy.integrate.solve_ivp.

Parameters:

  • _ – Time is inputted here by scipy.integrate.solve_ivp, but it is not used to calculate the Jacobian.

  • y – The current concentrations of each chemical.

Returns:

The Jacobian.

Return type:

np.ndarray

example

The simple chemical system:

\[A \xrightarrow{\text{k1}} B \xrightarrow{\text{k2}} C\]

can be modeled using the DRL class. First the reaction scheme should be written in a code friendly way:

reaction1 = ('k1', ['A'], ['B'])
reaction2 = ('k2', ['B'], ['C'])
reactions = [reaction1, reaction2]

Where the first element of each tuple is the name of the corresponding rate constant, the second element is a list containing all reactants, and the third element is a list containing all the products. If for example B split into C and byproduct D, we could write reaction2 as reaction2 = ('k2', ['B'], ['C', 'D'])

Lets assume that we know the rate constants belonging to this reaction.

rate_constants = {
    "k1": 0.2,
    "k2": 0.5,
}

We can create a prediction using the DRL.predict_concentration(). The class implements the method which determines the rate of change as a function of its current state, and a method which calculates the Jacobian based on its state. Because we do not want to model an entire DRL experiment, solve_ivp is used instead of DRL.predict_concentration(). Internally, this function also calls solve_ivp.

import numpy as np
from scipy.integrate import solve_ivp
from delayed_reactant_labeling.predict import DRL

time = np.linspace(0, 20, num=2000)  # desire predictions at these timestamps
k1, k2 = rate_constants['k1'], rate_constants['k2']
A0 = 1

drl = DRL(rate_constants=rate_constants, reactions=reactions, output_order=['A', 'B', 'C'], verbose=False)
result = solve_ivp(
    drl.calculate_step,
    t_span=[time[0], time[-1]],
    y0=[A0, 0, 0],
    method='Radau',
    t_eval=time,
    jac=drl.calculate_jac)

However, also algebraic solutions for this specific chemical problem exist.

\begin{eqnarray} [A]_t = [A]_0 \cdot e^{-k_1t} \end{eqnarray} \begin{eqnarray} [B]_t = \frac{k_1}{k_2-k_1}[A]_0(e^{-k_1t}-e^{-k_2t}) \end{eqnarray} \begin{eqnarray} [C]_t = [A]_0[1-e^{-k_1t}-\frac{k_1}{k_2-k_1}(e^{-k_1t}-e^{-k_2t})] \end{eqnarray}

We can compare the algebraic solution to the modelled prediction as follows.

import matplotlib.pyplot as plt
kinetic_A = A0 * np.exp(-k1 * time)
kinetic_B = k1 / (k2 - k1) * A0 * (np.exp(-k1 * time) - np.exp(-k2 * time))
kinetic_C = A0 * (1 - np.exp(-k1 * time) - k1 / (k2 - k1) * (np.exp(-k1 * time) - np.exp(-k2 * time)))

fig, ax = plt.subplots()
ax.plot(time, result.y[0] / A0, label='A')
ax.plot(time, result.y[1] / A0, label='B')
ax.plot(time, result.y[2] / A0, label='C')
ax.plot(time, kinetic_A, color='k', linestyle=':', label='algebraic')
ax.plot(time, kinetic_B, color='k', linestyle=':')
ax.plot(time, kinetic_C, color='k', linestyle=':')
ax.legend()
fig.show()
_images/predict_prediction.png

It is clear that the model fits the data very well. Besides, it’s much easier to implement these few lines of code, instead of doing the mathematics. Furthermore, implementing more difficult problems only requires the addition of a few lines here, whereas solving the problem in an exact manner might become impossible.

This current system can also be converted to a DRL experiment by adding a labeled reactant A at a certain timestamp. This labeled reactant reacts in a identical manner to A, except that the corresponding products are also labeled. It could be implemented as follows:

reaction3 = ('k1', ['A-labeled'], ['B-labeled'])
reaction4 = ('k2', ['B-labeled'], ['C-labeled'])
reactions.extend([reaction3, reaction4])

drl = DRL(reactions=reactions, rate_constants=rate_constants)
prediction = drl.predict_concentration(
    t_eval_pre=np.linspace(0, 2.5, 2500),
    t_eval_post=np.linspace(2.5, 20, 17500),
    initial_concentrations={"A": 1},
    labeled_concentration={"A-labeled": 0.8},
    dilution_factor=0.8
)
ax = prediction.plot('time')
ax.set_xlabel("time")
ax.set_ylabel("concentration")
ax.figure.show()
_images/predict_drl_prediction.png