chap19.ipynb

Kernel: Python 3 (ipykernel)

Printed and electronic copies of Modeling and Simulation in Python are available from No Starch Press and Bookshop.org and Amazon.

Case Studies Part 2

Modeling and Simulation in Python

Copyright 2021 Allen Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

# download modsim.py if necessary

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
    
download('https://raw.githubusercontent.com/AllenDowney/' +
         'ModSimPy/master/modsim.py')

# import functions from modsim

from modsim import *

This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. Click here to access the notebooks: https://allendowney.github.io/ModSimPy/.

This chapter presents case studies where you can apply the tools we have learned so far to the glucose-insulin minimal model, an electronic circuit, a thermal model of a wall, and the interaction of HIV and T cells.

The Glucose Minimal Model

In the previous chapter we implemented the glucose minimal model using given parameters, but I didn't say where those parameters came from.

In the repository for this book, you will find a notebook, glucose.ipynb, that shows how we can find the parameters that best fit the data. You can download it from https://github.com/AllenDowney/ModSimPy/raw/master/examples/glucose.ipynb or run it on Colab at https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/glucose.ipynb.

It uses a SciPy function called leastsq, which stands for "least squares"; so-named because it finds the parameters that minimize the sum of squared differences between the results of the model and the data.

You can think of leastsq as an optional tool for this book. We won't use it in the text itself, but it appears in a few of the case studies.

The Insulin Minimal Model

Along with the glucose minimal model, Berman et al. developed an insulin minimal model, in which the concentration of insulin, $I$, is governed by this differential equation:

$$\frac{dI}{dt} = -k I(t) + \gamma \left[ G(t) - G_T \right] t$$

where

• $k$ is a parameter that controls the rate of insulin disappearance independent of blood glucose.

• $G(t)$ is the measured concentration of blood glucose at time $t$.

• $G_T$ is the glucose threshold; when blood glucose is above this level, it triggers an increase in blood insulin.

• $\gamma$ is a parameter that controls the rate of increase (or decrease) in blood insulin when glucose is above (or below) $G_T$.

The initial condition is $I(0) = I_0$. As in the glucose minimal model, we treat this initial concentration as a free parameter; that is, we'll choose it to fit the data.

The parameters of this model can be used to estimate $\phi_1$ and $\phi_2$, which are quantities that "describe the sensitivity to glucose of the first and second phase pancreatic responsivity". These quantities are related to the parameters as follows:

$$\phi_1 = \frac{I_{max} - I_b}{k (G_0 - G_b)}$$$$\phi_2 = \gamma \times 10^4$$

where $I_{max}$ is the maximum measured insulin level, and $I_b$ and $G_b$ are the basal levels of insulin and glucose.

In the repository for this book, you will find a notebook, insulin.ipynb, that contains starter code for this case study. Use it to implement the insulin model, find the parameters that best fit the data, and estimate $\phi_1$ and $\phi_2$. You can download it from https://github.com/AllenDowney/ModSimPy/raw/master/examples/insulin.ipynb or run it on Colab at https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/insulin.ipynb.

Low-pass Filter

The following circuit diagram (from https://commons.wikimedia.org/wiki/File:1st_Order_Lowpass_Filter_RC.svg) shows a low-pass filter built with one resistor and one capacitor.

A filter is a circuit that takes a signal, $V_{in}$, as input and produces a signal, $V_{out}$, as output. In this context, a signal is a voltage that changes over time.

A filter is low-pass if it allows low-frequency signals to pass from $V_{in}$ to $V_{out}$ unchanged, but it reduces the amplitude of high-frequency signals.

By applying the laws of circuit analysis, we can derive a differential equation that describes the behavior of this system. By solving the differential equation, we can predict the effect of this circuit on any input signal.

Suppose we are given $V_{in}$ and $V_{out}$ at a particular instant in time. By Ohm's law, which is a simple model of the behavior of resistors, the instantaneous current through the resistor is:

$$I_R = (V_{in} - V_{out}) / R$$

where $R$ is resistance in ohms ($\Omega$).

Assuming that no current flows through the output of the circuit, Kirchhoff's current law implies that the current through the capacitor is:

$$I_C = I_R$$

According to a simple model of the behavior of capacitors, current through the capacitor causes a change in the voltage across the capacitor:

$$I_C = C \frac{d V_{out}}{dt}$$

where $C$ is capacitance in farads (F). Combining these equations yields a differential equation for $V_{out}$:

$$\frac{d V_{out}}{dt} = \frac{V_{in} - V_{out}}{R C}$$

In the repository for this book, you will find a notebook, filter.ipynb, which contains starter code for this case study. You can download it from https://github.com/AllenDowney/ModSimPy/raw/master/examples/filter.ipynb or run it on Colab at https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/filter.ipynb. Follow the instructions to simulate the low-pass filter for input signals like this:

$$V_{in}(t) = A \cos (2 \pi f t)$$

where $A$ is the amplitude of the input signal, say 5 V, and $f$ is the frequency of the signal in Hz.

Thermal Behavior of a Wall

This case study is based on a paper that models the thermal behavior of a brick wall with the goal of understanding the "performance gap between the expected energy use of buildings and their measured energy use".

The following figure shows the scenario and their model of the wall:

On the interior and exterior surfaces of the wall, they measure temperature and heat flux (rate of heat flow) over a period of three days. They model the wall using two thermal masses connected to the surfaces, and to each other, by thermal resistors.

The primary methodology of the paper is a statistical method for inferring the parameters of the system (two thermal masses and three thermal resistances).

The primary result is a comparison of two models: the one shown here with two thermal masses, and a simpler model with only one thermal mass. They find that the two-mass model is able to reproduce the measured fluxes substantially better.

For this case study we will implement their model and run it with the estimated parameters from the paper, and then use leastsq to see if we can find parameters that yield lower errors.

In the repository for this book, you will find a notebook, wall.ipynb with the code and results for this case study. You can download it from https://github.com/AllenDowney/ModSimPy/raw/master/examples/wall.ipynb or run it on Colab at https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/wall.ipynb.

The paper this case study is based on is Gori, Marincioni, Biddulph, Elwell, "Inferring the thermal resistance and effective thermal mass distribution of a wall from in situ measurements to characterise heat transfer at both the interior and exterior surfaces", Energy and Buildings, 2017, available from http://modsimpy.com/wall2.

The authors put their paper under a Creative Commons license and made their data available at http://modsimpy.com/wall. I thank them for their commitment to open, reproducible science, which made this case study possible.

HIV

During the initial phase of HIV infection, the concentration of the virus in the bloodstream typically increases quickly and then decreases. The most obvious explanation for the decline is an immune response that destroys the virus or controls its replication. However, at least in some patients, the decline occurs even without any detectable immune response.

In 1996 Andrew Phillips proposed another explanation for the decline in this paper: "Reduction of HIV Concentration During Acute Infection: Independence from a Specific Immune Response", available from https://people.math.gatech.edu/~weiss/uploads/5/8/6/1/58618765/phillips1996.pdf).

Phillips presents a system of differential equations that models the concentrations of the HIV virus and the CD4 cells it infects. The model does not include an immune response; nevertheless, it demonstrates behavior that is qualitatively similar to what is seen in patients during the first few weeks after infection.

His conclusion is that the observed decline in the concentration of HIV might not be caused by an immune response; it could be due to the dynamic interaction between HIV and the cells it infects.

In the repository for this book, you will find a notebook, hiv_model.ipynb, which you can use to implement Phillips's model and consider whether it does the work it is meant to do. You can download it from https://github.com/AllenDowney/ModSimPy/raw/master/examples/hiv_model.ipynb or run it on Colab at https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/hiv_model.ipynb.