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:


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))

# 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')


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!

  • $\begingroup$ If you're not too invested in odeint, the new api is recommended and has more solvers available docs.scipy.org/doc/scipy/reference/generated/… $\endgroup$
    – user9794
    Sep 2, 2023 at 13:16
  • $\begingroup$ It seems like you already know how to solve for different values of q and beta. It's not clear how you would keep H constant as it is the state variable you are solving for. $\endgroup$
    – user9794
    Sep 2, 2023 at 13:25
  • $\begingroup$ It's a good idea to be specific about what is not working, so sharing a direct comparison of what result you get and what result you should get would help. Quickly running the code you have and looking at the paper, I suspect you might have some issue with units. The variables and parameters must be in a consistent unit system. The paper uses seconds, hours and days, so perhaps something is off in some of your parameters. You are integrating with time unit of hours so all quantities should be given in units of hours or per hour. $\endgroup$
    – user9794
    Sep 3, 2023 at 9:12

1 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.

  • $\begingroup$ Thanks for the answer. As i mentioned that I am trying to reproduce the results of the paper that I mentioned, the issue is if you go through figure 2a and 2b, when I run my code I am not able to reproduce the same result. Even when I change the value of beta to 0.1, I am not able to reproduce figure 2c and 2d. So I guess my model either might be missing some parameter or maybe I coded it the wrong way $\endgroup$ Sep 2, 2023 at 18:35
  • $\begingroup$ I'm not going to reproduce the paper, so I'm basing this on what you've written. Your code as written only solves for 1 species, so if the paper has a system of interacting components, then your current code will not account for their interactions. You need to write the model as a system of equations and solve them simultaneously. $\endgroup$
    – user9794
    Sep 3, 2023 at 8:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.