Process

• Prepared by:

  • Joy Zhang

• Covered topics:

  • 0. Recap

  • 1. Understanding Process/Processes/processes

  • 2. Making a simple biokinetic model

  • 3. Advanced features

• Video demo:

  • Joy Zhang

To run tutorials in your browser, go to this Binder page.

You can also watch a video demo on YouTube (part 1, part 2) (subscriptions & likes appreciated!).

---

0. Recap

Check out this video here about qsdsan, make sure you understand the UML (unified modeling language) diagram.

from IPython.display import Image
Image(url='https://lucid.app/publicSegments/view/c8de361f-7292-47e3-8870-d6f668691728/image.png', width=800)

This tutorial will focus on the Process module in qsdsan. This module was developed to provide modeling capacity for dynamic conversions of waste streams in a SanUnit. With this module, simulations can describe not only a system’s mass and energy flows at steady states but also the changes of these flows over time.

import qsdsan as qs
print(f'This tutorial was made with qsdsan v{qs.__version__}.')

Back to top

1. Understanding Process/Processes/processes

Similar to the Component module, the Process module includes three main classes: Process, Processes, and CompiledProcesses.

The Process class is used to model individual dynamic processes with user-defined stoichiometry and reaction kinetics. The Processes class is simply a collection of multiple Process objects. Compilation of Processes yields CompiledProcesses, which provides additional attributes of the collection. It is often used to implement (pseudo-)mechanistic models with parellel biological and/or physiochemical processes in dynamic simulations of treatment systems.

Image(url='https://lucid.app/publicSegments/view/2c231fa2-6065-46b9-83af-a790ce38b6c0/image.png', width=600)

In this section, we will focus on understanding the structures of these classes by going over some existing processes in qsdsan. Next section will guide you through the creation of your own Processes subclass.

1.1. Processes vs. processes

You might remember SanUnit vs. sanunits. Processes and processes have similar distinctions.

# If you check
qs.Process
# or
# qs.Processes
# or
# qs.CompiledProcesses

This means Process/Processes/CompiledProcesses is a class in the _process module (i.e., the codes defining theses classes are in the _process.py file within the qsdsan folder).

# If you check
qs.processes

This shows that processes is actually a folder. (Why? See here.) It contains the process models (i.e., subclasses of Process or CompiledProcesses) that have been included in qsdsan.

# To see the list of objects that can be directly imported from this folder
qs.processes.__all__

The list above include common process models like ASM1, ASM2d, ADM1, which are implemented as subclasses of CompiledProcesses. You can also directly import utility functions that create CompiledComponents objects corresponding to the state variables in those process models (i.e., create_<ModelName>_cmps) among others.

Decay and KineticReaction are classes for simple processes whose state variables can be solved analytically.

# You can see other attributes of the `qs.processes` folder with the `dir` function
set(dir(qs.processes)) - set(qs.processes.__all__)

While names with double “_” on both sides (e.g., __path__) are built-in Python attributes, those with a single leading “_” (e.g., _asm1) are individual scripts included in the folder.

# For example, `_asm1.py` is a script in the `processes` folder.
qs.processes._asm1

For the rest of this section, we will use activated sludge model no. 1 (ASM1) as an example to learn about the core classes in the _process module.

Back to top

---

ASM1

ASM1 is a mathematical model for the activated sludge processes in wastewater treatment systems. It uses 13 components to characterize wastewater and activated sludge, including biomass, organic substrates (COD), organic and inorganic nitrogens (N), dissolved oxygen (DO), alkalinity etc.

ASM1 includes 8 transformation processes of the components (e.g., hydrolysis, biomass growth and decay). The stoichiometry of each process considers conservations of COD, N, and electric charges. The process rates are dependent on environmental conditions and state variables (i.e., relevant component concentrations).

# Before we get to the subsections, let's get ready by loading ASM1-related objects in qsdsan
from qsdsan.processes import create_asm1_cmps, ASM1
cmps = create_asm1_cmps()
cmps.show()

This means the create_asm1_cmps function returns a CompiledComponents object corresponding to the 13 components in ASM1 besides water.

Note: Classes in the Process module cannot be used independent of the Component module. Therefore, it is common practice to create corresponding CompiledComponents objects and set thermo before creating Process instances.

# By default, the thermo is set with this `CompiledComponents` object upon its creation.
# We can verify that by calling the `get_thermo` function
qs.get_thermo()

# Next we need to create an instance of the ASM1 model
# We can see that `ASM1` is a subclass of `CompiledProcesses`, so it can be used for demonstration
# of the capabilities of the `_process` module.
ASM1.__bases__

# Let's start by reading the documentation of the `ASM1` class
?ASM1

# Without getting into the details of ASM1, we will leave all parameters at their default values
asm1 = ASM1()
asm1.show()

This shows that asm1 is composed of the 8 transformation processes.

Back to top

---

1.2. Process

The Process class is developed to model individual dynamic processes. Dynamic means that the rate of reaction is finite (i.e., the reaction is not instantaneous) and often changes with the state of the system. Each of the eight transformation process in ASM1 is described by a Process object.

# Let's take the 0th process in `asm1` as an example to learn about `Process`:
# p1 = asm1.tuple[0]

# You can also access the process using the typical `getattr` method because
# individual processes in a `CompiledProcesses` object are stored as attributes
p1 = asm1.aero_growth_hetero
p1.show()

---

Tip

<Object>.<Attribute> is essentially identical to getattr(<Object>, '<AttributeName>').

---

p1 is an instance of Process describing the growth of heterotrophic biomass under aerobic conditions. The show method prints p1’s defining characteristics: stoichiometry, reference, rate equation, parameters, and dynamic parameters. All of them are stored as the object’s attributes.

# For example, we can retrieve information on the stoichiometry of this process
p1.stoichiometry

Unlike what’s printed by the show method, a dictionary of evaluated stoichiometric coefficients is returned. The coefficients indicate that to produce 1g COD’s worth of X_BH, approximately 1.49g COD’s worth of S_S, 0.08g N’s worth of S_NH, and 0.49g S_O are consumed. The values of these stoichiometric coefficients are dependent on the parameter values in the formulas. We’ll see how stoichiometry will change in response to different parameter values later.

Note: The stoichiometry must be defined and interpreted in consistency with the componenets’ units of measure. You can check the units by calling <Component>.measured_as.

Note that X_BH has a stoichiometric coefficient of 1, which means X_BH is a “reference component” (i.e., the process rate calculated by the rate equation is also X_BH’s rate of production/consumption).

# This information can also be accessed by calling the `ref_component` property.
p1.ref_component

# Another defining characteristics of a process is its rate equation, which is stored as a
# property of the `Process` object
p1.rate_equation

The process rate here is a “Monod-type” function of the concentrations (in \(g/m^3\)) of the four components (\(S_{NH}\), \(S_O\), \(S_S\), and \(X_{BH}\)) involved in this process.

How do I know this function takes concentrations instead of, let’s say, mass flowrates as input? Answer is – we can tell by looking at the units of the parameters, which are defined upon initiation of the Process object.

# If we access the private attribute
p1._rate_equation

This shows \(\mu_H\), \(K_{NH}\), \(K_{OH}\), and \(K_S\) are parameters in this rate equation. Although the information on their units are not stored with the Process object, we can see their definitions in the documentation of the ASM1 class.

# ?ASM1

Since \(\mu_H\) are defined in \(d^{-1}\), and \(K_{NH}\), \(K_{OH}\), and \(K_S\) in \(g/m^3\), only when \(S_{NH}\), \(S_O\), \(S_S\), and \(X_{BH}\) are in \(g/m^3\) can this rate equation output process rate in the correct unit (i.e., \(d^{-1}\)). This is defined by the mathematical formulation of ASM1. It is also common for biokinetic models to define state variables as component concentrations.

# Now that we've understood the required input to the rate equation, we can try evaluating the
# process rate. Let's try with all component concentrations equal to 1.
import numpy as np
state = np.ones(len(cmps)+1)
state

This is a numpy array with 15 elements. The first 14 elements are component concentrations in the exact same order as cmps, while the last element indicates volumetric flowrate \(Q\). Although \(Q\) is not used in process rate calculation, it is often an important state variable for the modeling of reactor-scale dynamics. Therefore, it is common for a lot of rate functions to ask for an array of length len(cmps)+1 as input, but it is ultimately determined by how the rate function is coded.

# Note that the `rate_equation` attribute only stores the formula.
# The evaluation of process rate is done through the `rate_function` attribute, which is rendered
# automatically from the rate equation.
p1.rate_function(state)

This result indicates that at the component concentrations and flowrate defined by the state array, this process happens at a rate of 0.152 \(d^{-1}\) (approximately), which translates to a heterotrophic biomass growth rate of 0.152 \(gCOD/m^{3}/d\).

Note: Both stoichiometry and rate_equation are properties of the Process object. Changing the formulas of stoichiometric coefficients or rate equation are prohibited after initiation of the object.

As you may remember from the ASM1 documentation and the asm1.show() printout above, parameters and/or dynamic parameters are also important attributes for a Process object because they represent intrinsic characteristics of the process and affect the evaluation of stoichiometry and process rate.

# For `asm1` and the individual processes within `asm1`, parameters are stored as a dictionary
p1.parameters
# asm1.parameters
# p1.parameters is asm1.parameters

You may notice p1 and asm1 share the dictionary of parameters, although not all parameters included in this dictionary are used in the stoichiometry or process rate calculations for p1. This is to make sure the parameter values are consistent across different processes in asm1.

# If you update parameter values for `p1`, the same parameters will be updated accordingly for
# `asm1` and any other processes in `asm1`.
p1.set_parameters(mu_H=6.0, K_S=8.0, Y_H=0.8)
asm1.parameters

# If you evaluate process rate or stoichiometry again with the same input,
# you should now expect different output.
# p1.stoichiometry
p1.rate_function(state)

In addition to the attributes above, conserved_for is another useful attribute, expecially for biochemical processes. It dictates the conservation principles followed by the stoichiometric coefficients of the the process.

p1.conserved_for

We can see p1 is subject to conservations of COD, N, and electric charge, as defined by ASM1. This information is useful for automatic derivation of unknown stoichiometric coefficients upon initiation of the Process object. If numerical values are provided for all stoichiometric coefficients when creating the Process object, this information can then be used by check_conservation to check if these values satisfy the conservation principles.

# No return indicates that all materials in `conserved_for` are conserved.
# Otherwise, a warning or a `RuntimeError` will be raised.
p1.check_conservation()

You may have noticed p1’s “dynamic parameters” is empty. This is because all of the parameters for the calculation of stoichiometric coefficients are static values, i.e., they don’t change during simulation. If a process’s stoichiometric coefficient is a function of time or other state variables but doesn’t have a nice and simple mathematical form, we can create a DynamicParameter object and add it to the process. We will learn how to do that in section 3.

Back to top

1.3. Processes and CompiledProcesses

It is more often than not you’re working with more than one dynamic process, which often makes it impossible to derive an analytical solution to the dynamics or the steady state of the system. This is when Processes and CompiledProcesses come in handy.

As is mentioned above, Processes is simply a collection of multiple Process objects. Compiling a Processes object yields a CompiledProcesses object. Note that CompiledProcesses is also a subclass of Processes, i.e., it has all of Processes’s attributes and more.

qs.CompiledProcesses.__bases__

We know that ASM1 is a subclass of CompiledProcesses, which makes asm1 an instance of CompiledProcesses.

# Let's verify that
isinstance(asm1, qs.CompiledProcesses)

So, is there a way to “decompile” a CompiledProcesses object? Yes!

asm1 = qs.Processes(asm1.tuple)
asm1.show()

# Let's see what difference "decompiling" made
asm1.__dict__

This shows that a Processes object only stores individual Process objects as its attributes.

asm1.compile()
asm1.__dict__

In comparison, a CompiledProcesses object has many more attributes that are convenient for the application of dynamic process models, although no additional information was provided for compile.

# For example, the stoichiometric coefficients of all processes are compiled into a table that
# is in consistent format as a Petersen matrix
asm1.stoichiometry

# Similarly for the rate equations
asm1.rate_equations

# More importantly, the `rate_function` attribute of a `CompiledProcesses` now outputs an array
# with each element corresponding orderly to the individual processes
asm1.rate_function(state)

With this, the evaluation of all process rates can be done with one command and the output data type is convenient for matrix operations.

To calculate the overall rate of production for a component, we need to first estimate its rates of production (positive means being produced, negative means being consumed) by each process and then sum them up. In mathematical forms, let’s say \(a_{ij}\) represents the stoichiometric coefficient of component \(j\) in process \(i\) (i.e., value on the \(i\)th row and \(j\)th column of the stoichiometry matrix). Denote \(\rho_i\) as process \(i\)’s reaction rate, now the overall production rate of component \(j\) should be

\[r_j = \sum_i{a_{ij}\cdot\rho_i}\]

In matrix notation, this calculation can be neatly described as

\[\mathbf{r} = \mathbf{A^T} \mathbf{\rho}\]

where \(\mathbf{A}\) is the stoichiometry matrix and \(\mathbf{\rho}\) is the array of process rates.

# This matrix operation is already streamlined for `CompiledProcesses` objects, you can see the
# mathematical form of the rates of production as a function of component concentrations.
asm1.production_rates

# To evaluate the rates of production for all components, all you need to
# do is to call the `production_rates_eval` method.
asm1.production_rates_eval(state)

With these convenient evaluation methods, asm1 can now be used in SanUnit to model system dynamics (e.g., reactor effluent concentrations over time). We will introduce how to integrate process models with SanUnit subclasses and how to enable dynamic simulations in a separate tutorial.

Back to top

2. Making a simple biokinetic model

Congratulations! You’ve understood the structure and the core capacities of the _process module. 🎉

In this section, we will get some hands-on experience with the _process module. For demonstration purpose, we will implement a simple biokinetic model with the classes we learned in the previous section. It’s even better if you have your own mathematical model of dynamic processes in mind.

Let’s get started with defining the mathematical model.

---

The model

Let’s consider a simplified activated sludge system with 5 components besides water: - particulate organic substrate (\(X_S\), measured as COD) - soluble organic substrate (\(S_S\), measured as COD) - dissolved oxygen (\(O_2\), measured as itself) - dissolved CO2 (\(CO_2\), measured as C) - biomass (\(X_B\), measured as COD)

This model uses 3 dynamic processes to describe the transformations among the components: 1. Hydrolysis (i.e., \(X_S\) transforming into \(S_S\)) 2. Growth of biomass (i.e., \(X_B\) utilizing \(S_S\) and \(O_2\) to grow) 3. Decay of biomass (i.e., \(X_B\) dying and transforming into \(X_S\))

To mathematically characterize these processes, we need to define the stoichiometry and the process rate equation of each process. Let’s say, for process 1 and 3:

\[X_S \Rightarrow S_S,~~\rho_1 = k_{hyd}X_S\]

\[X_B \Rightarrow X_S,~~\rho_3 = b\cdot X_B\]

where \(k_{hyd}\) is the hydrolysis rate constant, \(b\) is the biomass decay rate constant, both in \(d^{-1}\cdot m^3\cdot g^{-1}\).

For biomass growth, assuming the yield of biomass is \(y_B\) [gCOD \(X_B\)/gCOD \(S_S\)], then the stoichiometry can be defined as:

\[\frac{1}{y_B}S_S + a_1O_2 \Rightarrow X_B + a_2CO_2\]

where \(a_1\) and \(a_2\) need to be determined by conservations of COD and carbon. Similar to other activated sludge models, we can use a Monod-type function to represent biomass growth rate, i.e.,

\[\rho_2 = \mu_B\cdot \frac{S_S}{K_S+S_S}X_B\]

---

2.1. Defining the components

For convenience, we will find similar components that have been defined in the default set. Note that we will need to consider conservations of COD and C in the stoichiometry of process 2. Therefore, we should make sure the Component objects we pick have appropriate values for COD content (i_COD) and carbon content (i_C).

# Load the default set of components
cmps_all = qs.Components.load_default()

# Copy select components from the default set and give them new IDs
X_S = cmps_all.X_B_Subst.copy('X_S')
S_S = cmps_all.S_F.copy('S_S')
O2 = cmps_all.S_O2.copy('O2')
CO2 = cmps_all.S_CO3.copy('CO2')
X_B = cmps_all.X_OHO.copy('X_B')

# Compile them into a new set. Don't forget the water!
cmps_bkm = qs.Components([X_S, S_S, O2, CO2, X_B, cmps_all.H2O])
cmps_bkm.compile()
cmps_bkm.show()

Note: You might’ve noticed H2O isn’t an active component in any process but is always included in the CompiledComponents set. This is because water is assumed to be the bulk liquid and is needed for the convergence of total volumetric flows in simulation regardless of the process model.

# Now we can check if their `measured_as` attributes are correctly set
# X_S.measured_as == 'COD'
# CO2.measured_as == 'C'
O2.measured_as == None

# Then you can check relevant `i_` properties of the components
# For example, O2 should have a negative COD content, or more specifically -1 gCOD/gO2
O2.i_COD

# A final preparation step is to set thermo with the components
qs.set_thermo(cmps_bkm)

It is common practice to wrap up the creation of model-specific CompiledComponents objects in a function for convenient reuse. The create_<ModelName>_cmps functions we’ve seen in the previous section is a good example.

Back to top

2.2. Defining the processes

With the components ready, we can now define individual dynamic processes in the model as Process objects.

# As always, let's take a look at the documentation, especially the init signature.
?qs.Process

For process 1 and process 3, the stoichiometric coefficients are known and the rate equations take simple forms with just one parameter.

# To define `Process` objects for these two processes, we just need to input their
# stoichiometries and rate equations

# Process 1: hydrolysis
pc1 = qs.Process(ID='hydrolysis',
                reaction='X_S -> S_S',
                ref_component='X_S',
                rate_equation='k_hyd*X_S',
                conserved_for=(),
                parameters=('k_hyd',))

Note that the information on stoichiometry is reflected in the reaction argument. ref_component is set to be X_S because its stoichiometric coefficient is -1.

No conservation principle is used in the initiation of pc1 (i.e., conserved_for is an empty tuple). This is fine, because there’s no unknown stoichiometric coefficient that relies on conservation principles to solve. However, this doesn’t mean the process isn’t subject to any conservation principle. The 1 to 1 stoichiometry of this process reflects the conservation of COD, as both \(X_S\) and \(S_S\) are measured as COD.

# We can verify this with the `check_conservation` method
# To do that, `conserved_for` needs to be set first
pc1.conserved_for = ('COD',)
pc1.check_conservation()

No return means pc1 conserves COD.

Note:conserved_for is treated as an iterable in the check_conservation method. So, when there’s only one conservation principle, we should make sure the entire string is recognized as one element in iteration.

# For example, if we do this, we will get an error message because the function interprets it as
# checking conservations of 'C', 'O', and 'D', respectively
# pc1.conserved_for = ('COD')
# pc1.check_conservation()

Note that the initiation also asks for parameters. This is to help identify any symbolic parameter in the stoichiometry and the rate equation.

# Upon initiation, the parameters are stored as symbols. We still need to set values to them
# before we can evalute process rate.
pc1.show()

# At this point, the initiation of process 3 should be quite straightforward.
# Here shows an alternative way to input stoichiometry
pc3 = qs.Process('decay',
                reaction={'X_B':-1, 'X_S':1},
                ref_component='X_B',
                rate_equation='b*X_B',
                conserved_for=('COD',),
                parameters=('b',))

# But how do we initiate process 2, which has known stoichiometric coefficients?
pc2 = qs.Process('growth',
                reaction='[1/y_B]S_S + [?]O2 -> X_B + [?]CO2',
                ref_component='X_B',
                rate_equation='mu_B * S_S/(K_S+S_S) * X_B',
                conserved_for=('COD', 'C'),
                parameters=('y_B', 'mu_B', 'K_S'))
pc2.show()

Because there’re two unknown stoichiometric coefficients, we need two conservation principles (i.e., two equations) to determine their values. It is convenient to solve for stoichiometric coefficients marked as ? by specifying conserved_for in the initiation of Process objects.

Back to top

2.3. Compiling the processes

# Now the final step is to compile the individual processes into a biokinetic model
bkm = qs.Processes([pc1, pc2, pc3])
bkm.compile()
bkm.show()

# Parameters in the stoichiometry and rate equations across all processes are compiled into
# a shared dictionary.
# This also means, we should absolutely avoid using the same symbol for different parameters in
# multiple processes.
bkm.parameters

# After setting parameter values, the model will be ready
bkm.set_parameters(k_hyd=3.0, y_B=0.8, mu_B=4.0, K_S=9.0, b=0.4)
bkm.stoichiometry

bkm.production_rates

Back to top

2.4. Batch mode for creating processes

It’d be very tedious to create individual Process objects one by one and compile them, especially for models with a large number of parallel processes and/or components. Therefore, the _process module includes convenient classmethods to create processes in batch.

# Stoichiometry and rate equations are usually described in a table format
from qsdsan.utils import load_data
df_bkm = load_data('assets/_bkm.tsv', index_col=0)
df_bkm

This table contains the stoichiometry matrix and the array of rate equations. If the table is organized with process IDs as row indices and component IDs as column names, we can use it as direct input to the Processes.load_from_file classmethod to create processes in batch.

# The same amount of information is still required to create the model.
bkm_batch = qs.Processes.load_from_file(
    path='assets/_bkm.tsv',
    conserved_for=('COD', 'C'),
    parameters=('k_hyd', 'y_B', 'mu_B', 'K_S', 'b'),
    compile=True
)
bkm_batch.show()

# You can see `bkm_batch` is equivalent to the `bkm` we created above
bkm_batch.stoichiometry
# bkm_batch.rate_equations
# bkm_batch.parameters

# The reference component of each process is inferred from its stoichiometry
bkm_batch.decay.ref_component

# `conserved_for` now applies to all processes,
# the following will trigger an error
# bkm_batch.decay.check_conservation()

This means the decay process defined above conserves COD but not carbon, suggesting that we could revisit the i_C values for \(X_B\) and \(X_S\) or consider whether there are other carbon-containing products in this process.

# `i_C` values for each component
dict(zip(cmps_bkm.IDs, cmps_bkm.i_C))

Back to top

3. Advanced features

qsdsan aims to provide a flexible and user-friendly platform for process modeling. Therefore, many advanced features are included in the _process module to streamline the implementation of process models with various mathematical complexity.

Let’s again use the model above as an example to learn about these features.

3.1. DynamicParameter

What if the yield of biomass growth is dependent on the abundance of soluble organic substrate? Assume this relation can be described as

\[y_B = y_{max}\cdot\sqrt{\frac{S_S}{X_S + S_S}}\]

In this case, we will need to define \(y_B\) as a DynamicParameter object because its value changes with system state during simulation.

@pc2.dynamic_parameter(symbol='y_B', params={'y_max':0.8})
def yB_eval(state_arr, params):
    X_S, S_S = state_arr[:2]  # assuming the first 6 elements are component concentrations
    y_max = params['y_max']
    y_B = y_max * (S_S/(X_S+S_S))**(1/2)
    return y_B

The function yB_eval must take no more than two positional arguments, with the first one expecting an array of state variable values and the second one expecting a dictionary of additional inputs required for the evaluation of this dynamic parameter.

Previous tutorials have covered the use of decorator. You can revisit this tutorial for more details.

# At this point, the value for `y_B` is not updated yet, since we haven't evalutated it
# with input of component concentrations
bkm.parameters

# But the list of dynamic parameters have been updated
pc2.show()

# `y_B` is now a `DynamicParameter` object stored in the `_dyn_params` attribute of the process
type(pc2._dyn_params['y_B'])

# Assuming component concentrations are all 1.
state_bkm = np.ones(6)

# This function evaluates the dynamic parameters
bkm.params_eval(state_bkm)
bkm.parameters

# Then the stoichiometry is updated accordingly
bkm.stoichiometry

Note that in simulations, you don’t need to call the params_eval method explicitly unless you would like to capture the dynamics of parameter values. The evaluation of dynamic parameters is performed within production_rates_eval (i.e., the calculation of overall production rates of the components).

Back to top

3.2. Kinetics and MultiKinetics

Kinetics and MultiKinetics are classes that allow defining process rates as functions.

Kinetics is a subclass of DynamicParameter, because process rate in essence is also a function of state variables and additional static parameters.

# The use of `Kinetics` to define process rate is also similar
# For example, for the decay process
@pc3.kinetics(parameters={'b':0.4})
def rho3_eval(state_arr, params):
    X_B = state_arr[4]
    b = params['b']
    return b*X_B

pc3.rate_function

This may seem trivial, but it’s helpful when the form of process rate gets more complex. For example, the rate of substrate uptake in anaerobic digestion model no.1 is dependent on pH, which is dynamically determined by a linear equation system and thus cannot be expressed neatly as a function of state variables.

MultiKinetics serves a similar purpose but calculates the rate of multiple processes at onces. It is designed to be used with CompiledProcesses objects. This is convenient for models with similar process rate algorithms across multiple processes in the sense that vectorizing the shared algorithm can often improve computational speed in simulation. Another application of MultiKinetics is models sharing intermediate variables (e.g., pH) across multiple process rate calculation.

Let’s assume the reaction rates of all three processes in the model above are now dependent on the abundance of oxygen. Here are the new rate equations:

\[\rho_1 = k_{hyd}X_S\cdot \sqrt{\frac{O_2}{O_2+K_{O2}}}\]

\[\rho_2 = \mu_BX_B\cdot \frac{S_S}{K_S+S_S} \sqrt{\frac{O_2}{O_2+K_{O2}}}\]

\[\rho_3 = b\cdot X_B\cdot (1.5- \sqrt{\frac{O_2}{O_2+K_{O2}}})\]

How do we use MultiKinetics to implement this change?

# Define the function
def rhos_eval(state_arr, params):
    X_S, S_S, O2, CO2, X_B = state_arr[:5]
    I_O2 = (O2/(O2+params['K_O2']))**(1/2)
    rho1 = params['k_hyd']*X_S*I_O2
    rho2 = params['mu_B']*X_B*I_O2
    rho3 = params['b']*X_B*(1.5-I_O2)
    return np.array([rho1, rho2, rho3])

# Set it as the rate function of the `CompiledProcesses` object
bkm.set_rate_function(rhos_eval)

# Assign parameter values
bkm.rate_function.set_params(K_O2=2.0, **bkm.parameters)

# The evaluation of the rate function is the same as before
bkm.rate_function(state_bkm)

Back to top