Grammar

Grammar#

gotranx uses a domain specific language (DSL), that is implemented using lark, and you can find the complete grammar here.

In this section, we will try to describe the grammar, and how you can write your own .ode files to define your ODEs using this grammar.

Basic ODE file#

An example .ode file called lorentz.ode could look as follows

# Define parameters
parameters(sigma=12.0, rho=21.0, beta=2.4)

# Define states
states(x=1.0, y=2.0, z=3.05)

# Define expressions
dx_dt = sigma * (y - x)
a = rho - z
dy_dt = x * a - y
dz_dt = x * y - beta * z

We will now go through the basic building block for creating an ODE

parameters allows users to define parameters in the ODE. These are constant that the a user can provide as arguments. For example
```
parameters(sigma=12.0, rho=21.0, beta=2.4)
```
will create three parameters sigma, rho and beta with a default value of 12, 21.0 and 2.4 respectively.
states defines the state variables which are the variables you want to solve for. For example
```
states(x=1.0, y=2.0, z=3.05)
```
will define three state variables x, y and z who’s default initial condition will be set to 1, 2 and 3.0 respectively.
expressions allows users to define mathematical expressions using parameters and other variables. You can also define new variables as an expression. In this section you also need to specify the the state derivatives of all the state variables. For the state derivatives you need the variables name to be d{STATE_NAME}_dt, e.g dx_dt for the state derivative of x. In this case we could for example have
```
dx_dt = sigma * (y - x)
a = rho - z
dy_dt = x * a - y
dz_dt = x * y - beta * z
```
In this context, a will be referred to as an intermediate variable.
Comments can be added on separate lines by starting the line #, for example
```
# Define parameters
```

Variables can only be defined once#

An important thing to note is that once you have defined a variable you cannot redefine it, so for example the following expression is not valid

x = 1
y = 2 * x
x = 3

Operators and symbols#

Mathematical operators#

The following unary mathematical operators are supported

exp
cos
sin
tan
acos
asin
atan
abs
floor
ln
log
sqrt

The following binary mathematical operators are supported

Mod

For example

x = abs(1.0 - exp(4))

Numbers#

A basic number can be written as e.g

x = 1             # Will be casted to an integer
y = 1.0           # Will be casted to a float
z = 1 / 4         # Will be cased to a rational number
w = (2 * 3) / 3   # Will be cased to an expression

It is also possible to use scientific notation using either uppercase of lowercase e, e.g the following numbers are equivalent

x = 100
x = 1e2
x = 1E2

Similarly you can also use negative exponents, e.g the following numbers are equivalent

x = 0.01
x = 1e-2
x = 1E-2

Mathematical constants#

There is currently only one reserved constant with is pi. This means that you can e.g write

x = cos(2 * pi)

without the need to define pi in advance

Conditional operators and boolean operators#

The following boolean operators are supported

Lt which is equivalent to <
Gt which is equivalent to >
Le which is equivalent to <=
Ge which is equivalent to >=
Eq which is equivalent to ==
Not
And
Or

It is also possible to use the conditionals using the Conditional function (there is also a continuous version called ContinuousConditional), which has the following syntax

Conditional(boolean_expression, true_value, false_value)

For example the Heaviside step function

\[\begin{split} H(x) = \begin{cases} 1, \;\; x \geq 0, \\ 0, \;\; x < 0, \end{cases} \end{split}\]

can be defined as follows

H = Conditional(Geq(x, 0), 1, 0)

This would be similar to implemented an if-else statement, e.g

if x >= 0:
  H = 1
else:
  H = 0

We could also have a conditional with three value, e.g

\[\begin{split} H(x) = \begin{cases} 1, \;\; x > 0, \\ 0.5 \;\; x = 0, \\ 0, \;\; x < 0, \end{cases} \end{split}\]

which would be equivalent to and if-elif-else statement, e.g

if x > 0:
  H = 1
elif x == 0:
  H = 0.5
else:
  H = 0

but in this case we would need to nest conditionals together, for example

H = Conditional(Gt(x, 0), 1, Conditional(Eq(x, 0.0), 0.5, 0.0))

or another solution could be

H = Conditional(Gt(x, 0), 1, Conditional(Lt(x, 0.0), 0.0, 0.5))

Implementing `min` and `max`#

The operations min and max can be implemented in terms of Conditional, because

f = min(x, y)

is equivalent to

if x <= y:
  f = x
else:
  f = y

So f = min(x, y) could be implemented as

f = Conditional(Le(x, y), x, y)

Similarly f = max(x, y) can be implemented as

f = Conditional(Ge(x, y), x, y)

Advanced usage#

Adding components, units and descriptions#

If you have a large system of ODE’s it might be good to group parts that together in the same group, which we will refer to as components.

states("membrane",
V=ScalarParam(-87, unit="mV", description="")
)

states("sodium_channel_h_gate",
h=ScalarParam(0.8, unit="1", description="")
)

states("sodium_channel_m_gate",
m=ScalarParam(0.01, unit="1", description="")
)

states("potassium_channel_n_gate",
n=ScalarParam(0.01, unit="1", description="")
)

parameters("membrane",
Cm=ScalarParam(12.0, unit="uF", description="")
)

parameters("leakage_current",
E_L=ScalarParam(-60.0, unit="mV", description=""),
g_L=ScalarParam(75.0, unit="uS", description="")
)

parameters("sodium_channel",
E_Na=ScalarParam(40.0, unit="mV", description=""),
g_Na_max=ScalarParam(400000.0, unit="uS", description="")
)

expressions("sodium_channel_h_gate")
alpha_h = 170.0*exp((-V - 1*90.0)/20.0) # S/F
beta_h = 1000.0/(exp((-V - 1*42.0)/10.0) + 1.0) # S/F
dh_dt = alpha_h*(1.0 - h) - beta_h*h

expressions("sodium_channel_m_gate")
alpha_m = (100.0*(-V - 1*48.0))/(exp((-V - 1*48.0)/15.0) - 1*1.0) # S/F
beta_m = (120.0*(V + 8.0))/(exp((V + 8.0)/5.0) - 1*1.0) # S/F
dm_dt = alpha_m*(1.0 - m) - beta_m*m

expressions("potassium_channel_n_gate")
alpha_n = (0.1*(-V - 1*50.0))/(exp((-V - 1*50.0)/10.0) - 1*1.0) # S/F
beta_n = 2.0*exp((-V - 1*90.0)/80.0) # S/F
dn_dt = alpha_n*(1.0 - n) - beta_n*n

expressions("potassium_channel")
g_K1 = 1200.0*exp((-V - 1*90.0)/50.0) + 15.0*exp((V + 90.0)/60.0) # uS
g_K2 = 1200.0*n**4.0 # uS
i_K = (V + 100.0)*(g_K1 + g_K2) # nA

expressions("sodium_channel")
g_Na = g_Na_max*(h*m**3.0) # uS
i_Na = (-E_Na + V)*(g_Na + 140.0) # nA

expressions("leakage_current")
i_Leak = g_L*(-E_L + V) # nA

expressions("membrane")
dV_dt = (-(i_Leak + (i_K + i_Na)))/Cm # mV