PhysioModeler

PhysioModeler is designed as an educational tool to help model and simulate physiological systems¹, as represented by ordinary differential equations (ODEs). The design goal for this package is an explicit, easy to use, easy to read interface.²

A Model is defined by at least two pieces of information: what are the state components, and how do they change? A state component generally represents a physical quantity, like volume or distance. The rate of change of the state are described by differential equations.

pendulum_model = Model(
    dynamics=...,          # describes how the state change
    state_components=...,  # which state components to track
    parameters=...,        # model parameters
    inputs=...,            # inputs to the model
)

Example: undamped pendulum

A classic example is the humble undamped pendulum. If you know the position and angular velocity (i.e. how fast the pendulum is swinging) at any point in time, you can predict the movement of the pendulum in the future. Therefore, we track the angle \(\theta\) and angular velocity \(\frac{d\theta}{dt}\) as state.

The way a pendulum swings can be modelled using the second derivative of the angle the pendulum makes with a vertical line, known as the angular acceleration \(\frac{d^2\theta}{dt}\).

\[ \frac{d^2\theta}{dt^2} = -\frac{g}{l}\sin\theta \]

The angular acceleration \(\frac{d^2\theta}{dt^2}\) describes how the angular velocity \(\frac{d\theta}{dt}\) changes over time, while the angular velocity describes how the angle \(\theta\) changes over time. The Python-function below receives the state components \(\theta\) ('theta') and \(\frac{d\theta}{dt}\) ('dtheta'), and returns how they change over time, respectively as \(\frac{d\theta}{dt}\) ('dtheta') and \(\frac{d^2\theta}{dt^2}\) ('ddtheta').

import math

def update_angle_angular_velocity(time, state, inputs, parameters):
    l = parameters['length']
    g = inputs['gravity']

    return {
        "dtheta": state['dtheta'],
        "ddtheta": -g / l * math.sin(state['theta'])
    }

We now have the minimal amount of information requires to create a Model. However, there are two pieces of information missing: the parameter \(l\) for the length of the pendulum, and the gravity \(g\) that pulls on the pendulum. These can be provided as parameters and inputs. Parameters are constant values that describe characteristics of the model itself, while inputs describe external influences which might change over time. In this case, we assume gravity is constant.

pendulum_model = Model(
    state_components=["theta", "dtheta"],        # which state components to track
    function=update_angle_angular_velocity,  # how the state change
    parameters={"l": 1},                     # model parameters
    inputs={"g": 9.81},                      # external values
)

We can now run a simulation of this model. If we don't tell it otherwise, the simulation will start with \(\theta=0\,\mathrm{rad}\) and \(\frac{d\theta}{dt}=0\,\mathrm{rad}/\mathrm{s}\). Since a static pendulum is not much fun at all, let's have it start at an angle of \(1\,\mathrm{rad}\).

simulation_result = pendulum_model.run_simulation(
    time=60, 
    initial_state={"theta": 1}
)

The simulation result is a pandas DataFrame, which supports easy plotting of the results.

simulation_result.plot(subplots=["theta", "dtheta", "ddtheta"])

Including the parameters and inputs in the model definition (instead of depending on external dictionaries or hard-coded values) allow for easy comparison of different scenarios:

shorter_string_result = pendulum_model.run_simulation(
    time=60, 
    initial_state={"theta": 1},
    parameters={"l": 0.5}
)
moon_result = pendulum_model.run_simulation(
    time=60, 
    initial_state={"theta": 1},
    inputs={"g": 1.625}
)

---

1. Although we designed this software for physiological systems, it should be able to cope with a wide range of ODE models. ↩

2. We mostly try to adhere to lines 2 and 7 of the Zen of Python: "Explicit is better than implicit", and "Readability counts". ↩