Important: The unit system is a feature under development. The following section might not work properly yet.
Here, we show some examples using the unit system in ecell4. This feature requires Python library, pint. Install pint before running this example as follows: pip install pint. See also the website https://pint.readthedocs.io/en/latest/.
In [1]:
%matplotlib inline
from ecell4 import *
In [2]:
with species_attributes():
A | B | C | {'D': 1, 'radius': 0.005}
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
The species_attributes section defines a diffusion constant and radius of Species, A, B and C. For example, the diffusion rate of A is 1, and its dimensionality is expected to be [length**2/time]. However, what is the scale? Is it meter? Or mile?
Once the base units are determined, e.g. micrometer as [length] and second as [time], all the units must be consistent within the model. The second order rate consant must have dimensionality [1/(substance/length**3)/time], which is micrometer**3/item/second. Thus, when the parameter is given as 1/molar/min in some literature, you have to translate it by yourself.
In [3]:
show(m)
In [4]:
from ecell4.extra.unit import getUnitRegistry
In [5]:
ureg = getUnitRegistry()
Q_ = ureg.Quantity
with species_attributes():
A | B | C | {'D': Q_(1, 'um**2/s'), 'radius': Q_(0.005, 'um')}
with reaction_rules():
A + B == C | (Q_(0.01, '1/(item/um**3)/s'), Q_(0.3, '1/s'))
m = get_model()
First, create your own unit system,ureg, by using ecell4.extra.unit.getUnitRegistry. With this UnitRegistry, you can make a quantity with its unit as ureg.Quantity(value, unit). (Please be careful about the type of Quantity. It looks same with Quantity given by pint, but is slightly changed in ecell4 though all the original functionality in pint is availble even in ecell4. Please not use ureg = pint.UnitRegistry().)
In [6]:
help(ureg.Quantity)
The default base units are meter for [length], second for [time], and item (which means the number of molecules) for [substance]. When you change the default base unit, do like ureg = getUnitRegistry(length='micrometer').
In [7]:
show(m)
Now you can provide quantities in any unit regardless of the base units.
In [8]:
ureg = getUnitRegistry()
Q_ = ureg.Quantity
with species_attributes():
A | B | C | {'D': Q_(1e-8, 'cm**2/s'), 'radius': Q_(5, 'nm')}
with reaction_rules():
A + B == C | (Q_(6.02214129, '1/uM/s'), Q_(18, '1/min'))
m = get_model()
show(m)
You can operate quantities, and make a new quantity. See https://pint.readthedocs.io/en/latest/ for more details.
In [9]:
volume = Q_(1, 'fL')
conc = Q_(100, 'nM')
print((volume * conc).to('item'))
In addition to the model creation, run_simulation (and ensemble_simulations) also supports the unit system.
In [10]:
run_simulation(Q_(0.1, 'min'), y0={'C': Q_(60, 'item')}, volume=Q_(1, 'fL'), model=m, solver='ode')
Even if you change the base units, the behavior of simulations is kept consistent. In the following example, base units are rescaled to micrometer and minute with no change in the modeling section.
In [11]:
ureg = getUnitRegistry(length='micrometer', time='minute')
Q_ = ureg.Quantity
with species_attributes():
A | B | C | {'D': Q_(1e-8, 'cm**2/s'), 'radius': Q_(5, 'nm')}
with reaction_rules():
A + B == C | (Q_(6.02214129, '1/uM/s'), Q_(18, '1/min'))
m = get_model()
show(m)
run_simulation(Q_(0.1, 'min'), y0={'C': Q_(60, 'item')}, volume=Q_(1, 'fL'), model=m, solver='ode')
As default, the dimensionality of units is checked. When you turn off unit checking, call ecell4.extra.unit.use_strict(False) before using those features.
For example, the dimensionality of a diffusion constant must be [length**2/time]. When you give a unit with a wrong dimensionality, an exception would be thrown like:
In [12]:
try:
with species_attributes():
A | {'D': Q_(1.0, 'um/s')}
except ValueError as e:
print('{}: {}'.format(e.__class__.__name__, e))
Similarly, a kinetic rate constant of reactions is verified based on the order of the reaction. The first order reaction rate should have [1/time], and the second order should have [l/(substance/length**3)/time] in volume.
try: with reactionrules(): A + B > C | Q(0.3, '1/s') except ValueError as e: print('{}: {}'.format(e.class.name, e))
Additionally, rate law representations accept quantities with a unit too. See the example below:
In [13]:
ureg = getUnitRegistry()
Q_ = ureg.Quantity
with reaction_rules():
S > P | Q_(1.0, 'uM/s') * S / (Q_(100, 'nM') + S)
m = get_model()
show(m)
Here, the reaction above has two quantities, Vmax = Q_(1.0, 'uM/s') and Km = Q_(100, 'nM'). First, Km must have the same dimensionality with S, which is [concentration].
In [14]:
try:
with reaction_rules():
S > P | Q_(1.0, 'uM/s') * S / (Q_(100, 'nM/s') + S)
except ValueError as e:
print('{}: {}'.format(e.__class__.__name__, e))
Secondly, the dimensionality of a rate equation must be [concentration/time]. Therefore, the dimensionality of Vmax should be [concentration/time] too.
In [15]:
try:
with reaction_rules():
S > P | Q_(1.0, '1/s') * S / (Q_(100, 'nM') + S)
except RuntimeError as e:
print('{}: {}'.format(e.__class__.__name__, e))
When you give a value with no unit, it is regarded as dimensionless.
In [16]:
with reaction_rules():
S > P | 10.0 * Q_(0.1, 'uM/s') * S**2 / (Q_(100, 'nM')**2 + S**2)
m = get_model()
show(m)