Modeling Homeostatic Regulation with Python: How to Simulate Dynamic Systems

Question

Hello Computational Science community,

I am working on modeling a dynamic system that involves homeostatic regulation in Python and I could use some guidance on how to implement it. Here's a brief overview of the problem:

Background:

I'm trying to simulate a biological system where different metabolites (chemical species) are produced and cleared from the system, and their concentrations are regulated by homeostatic mechanisms. Each metabolite has its own production rate ($q$), and they can interact with each other through obstruction effects (beta) and circadian modulation. Basically I want to reproduce the results of the following paper: "Cortical waste clearance in normal and restricted sleep with potential runaway tau buildup in Alzheimer’s disease"

Model Equations:

The key equations governing this system are as follows:

Clearance Rate ($\chi_{\sigma_i}(t)$):

$$ \chi_{\sigma_i}(t) = \frac{\chi_{{\sigma_i}_0} }{ (1 + \sum_{j=1}^n \beta_{\sigma_j} H_j(t))} \chi_{\sigma_i}(t) $$ represents the clearance rate of species i at time $t$. $\chi_{{\sigma_i}_0}$ is the initial clearance rate of species i. beta_sigma_j is the obstruction rate due to species j. $H_j$ represents the homeostatic pressures of species j.

Homeostatic Pressure ($H_i(t)$): $$ \frac{dH_i(t)}{dt} = q_{\sigma_i} - \chi_{\sigma_i}(t) H_i(t) $$ represents the homeostatic pressure of species i at time t. $q_{\sigma_i}$ is the production rate of species i. $\chi_{\sigma_i}(t)$ is the clearance rate of species i at time t.

sigma in all the cases can take two forms: either sleep(s) or wake(w). so the input for sleep or wake has to be given by the user in the model.

Python Implementation: I'm trying to implement this system in Python using numerical methods. I have already set up the differential equations and used the odeint function from SciPy to simulate the system for a single species with varying production rates and circadian modulation. Here's my code so far:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint

# Metabolite parameters
q_values = np.array([7e-6, 65e-6, 140e-6])  # Production rates for q_sigma1
beta_value = 0  # Obstruction effect
chi_0 = 0.06

# Define the differential equations
def model(H, t, q, beta):
    chi = chi_0 / (1 + np.sum(beta * H))
    dH_dt = q - chi * H
    return dH_dt

# Simulation parameters
total_time = 24 * 8  # One week in hours
time_points = np.linspace(0, total_time, 1000)

# Initialize arrays to store results
H_initial = 0
results = []

# Simulate for different production rates
for q in q_values:
    result = odeint(model, H_initial, time_points, args=(q, beta_value))
    results.append(result)

# Plot the results
plt.figure(figsize=(10, 6))
for j, q in enumerate(q_values):
    plt.plot(time_points, results[j][:, 0], label=f'q = {q:.2e}')

plt.xlabel('Time (hours)')
plt.ylabel('Homeostatic Drive')
plt.title('Simulation for Single Metabolite with Varying Production Rate and Circadian   Modulation')
plt.legend()
plt.show()`

Question:

Now, I would like to generalize this model to simulate it for multiple species (n > 1) and different combinations of beta and q values while keeping other parameters constant (for example H=0 or H=10 and beta=0, or beta=0.1, or beta=0.2 etc). Additionally, I want to add the circadian drive ($C = \sin(w t)$) to the $\chi_{\sigma_i}(t)$ equation as what the authors have done in the above paper.

Could someone please guide me on how to extend my current code to achieve this? I'm particularly interested in how to set up the system for multiple species and introduce the circadian modulation.

Any help or suggestions would be greatly appreciated. Thank you!

Answer 1

To solve for multiple species you just need to make your model function return a vector where each entry is the derivative of one variable. The example in the documentation shows this for a system of two equations https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint

def model(y, t):
    # y is now a vector of state variables, y[0], y[1] and so on
    dy0_dt = ... # derivative of y[0]
    dy1_dt = ... # derivative of y[1]
    return [dy0_dt, dy1_dt]

It seems like you could just have two model functions, one for sleep and one for wake, and just use the appropriate one when you need. If you know explicitly when sleep and wake periods start and end, you can just switch equations and solve each interval with the appropriate model using the final value of the previous interval as the initial value when you solve the next interval.