Recode: Extraordinary intelligence and the care of infants

We recode the model of the article "Extraordinary intelligence and the care of infants" (10.1073/pnas.1506752113) by Steve Piantadosi and Celeste Kidd. The pdf is available here. Here, we only succinctly describe the model. You should consult the original article for details and for the rationale behind the model's choices.

The spirit of this notebook is to use simple code that is easy to understand and modify. This notebook requires no specific knowledge beyond a basic grasp of the Python language. We show all the code of the model, without relying on any library beyond numpy. Only the plotting, using the bokeh library, have been abstracted away in the graphs.py file. We employ the reproducible library to keep track of the computational environment and foster reproducibility.

A citable version of this notebook is available at figshare. You can contact me for questions or remarks at fabien@benureau.com.

In [1]:
import numpy as np
import graphs

np.random.seed(0)                                                                               # repeatable results.

The reproducible library fosters reproducibility. It collects information about the OS, Python, CPU and the git repository if present, so it can be displayed in the last cell of this notebook.

In [2]:
import reproducible
context = reproducible.Context()                          # the Context instance collects data about the environment.
try:
    context.add_repo(path='.', allow_dirty=True, diff=True)             # collect the hash of the current git commit.
except reproducible.RepositoryNotFound:
    pass

Model Equations

Piantosi and Kidd's model ties together three main parameters:

  • $R$, the size of the adult brain.
  • $T$, the duration of the gestation period.
  • $I_p$, a quantification of the intelligence of the parents.

The size of the child's head at age $t$, $g(t, R)$, follows a Gompertz growth curve, with $b$ and $c$ free parameters, fitted to 37.1 and 0.42 respectively.

$$g(t, R) = Re^{-be^{-ct}}$$
In [3]:
def g(t, R, b=37.1, c=0.42):
    """Size of the head at time t, given an adult size R."""
    return R * np.exp(-b * np.exp(-c * t)) # if you modify this function, you must modify the solve_MT function below.

Because a large head means a more difficult and dangerous birth, the probability to survive childbirth decreases (sigmoidally, here the function ϕ) when the head size at birth exceed a fixed parameter $V$, fitted to 5.48 cm. $T$ is the duration of the gestation period.

$$P(\textrm{survive childbirth}\,|\,T,R) = \phi(V - g(T,R))$$
In [4]:
def ϕ(z):
    """Sigmoid function"""
    return 1/(1 + np.exp(-z))

def P_born(T, R, V=5.48, g=g):
    """Probability to survive birth"""
    return ϕ(V - g(T, R))

The probability to survive adulthood is tied to the time after birth to reach maturity $M$ ($M$ solves $g(M + T, R) = 0.99R$), and the intelligence of the parents $I_p$.

$$P(\textrm{survive to adulthood}\,|\,M,I_p) = e^{-M(\gamma/I_p)}$$

Here, the free parameter $\gamma$, capturing the rate of mortality, is fitted to 0.4.

In [5]:
def P_adult(M, I_p, γ=0.4):
    """Probability to survive to adulthood"""
    return np.exp(-max(0, M) * γ / I_p)

The article assumes that $I_p$ is equal to the brain radius $R$ [1] for Figure 1 and 2A. Figure S2 of the supplementary material of the article explores other possible choices (area, volume, log volume). To try other relationships, modify the following function:

In [6]:
def I(R):
    """Return the intelligence (of the parents) as a function of R"""
    return R

Figure 1: Child Growth, Birth and Childhood Survival

We reproduce Figure 1. The continuous line correspond to $R$ = 8.4 cm, and the dashed line to $R$ = 4.2 cm. You can modify the latter if you are running the jupyter notebook version (the slider will not appear in the html one) by using the slider bellow.

In [7]:
ts = np.linspace(0, 25, 251)                                     ## i.e., ts = [0.0, 0.1, 0.2, ..., 24.8, 24.9, 25.0] 
fig1_data = graphs.fig1(ts, g, P_born, P_adult, I, R=8.4/2)
In [8]:
def update_fig1(R):
    graphs.update_fig1(fig1_data, g, P_born, P_adult, R)
graphs.interact(update_fig1, R=graphs.FloatSlider(min=0.1,max=20.0,step=0.1,value=4.2))

If the slider has no effect, rerun the last two code cells.

Figure 2A: Fitness Landscape

To reproduce Figure 2A, we need to compute $M$, which solves $g(M + T, R) = 0.99R$, the solution of which does not actually depends on the value of $R$. To speed up the computation, we use an analytical solution, only valid for $g(t, R) = Re^{-be^{-ct}}$.

To allow for arbitrary modification of the $g$ function, we also provide a simple dichotomy method to solve $M$.

In [9]:
def solve_MT(R, g, b, c):
    """Return M+T, with M+T solving g(M+T, R) == 0.99*R."""
    return -np.log(-np.log(0.99) / b) / c       # closed-form solution. Not valid if you modify the g function above.
                                                # In this case, comment this line, the code below is general-purpose.
    low, up = 1e-3, 25
    while up-low > 1e-3:                                                                # simple dichotomy algorithm.
        middle = 0.5*(up + low)
        if g(middle, R, b=b, c=c) < 0.99*R:
            low = middle
        else:
            up = middle
    return 0.5*(up+low)
In [10]:
K = 400                                                                              # K resolution of the landscape.
Ts = np.linspace(  0, 30, K+1)                                              # birth age, K+1 points between 0 and 30. 
Rs = np.linspace(0.1, 10, K)                                  # brain size (radius, cm), K points between 0.1 and 10.

def probability_matrix(Ts, Rs, γ=0.4, V=5.48, b=37.1, c=0.42):
    """Return the matrix of the probabilities to survive until adulthood."""
    D = []
    for R in Rs:
        D.append([])
        MT = solve_MT(R, g, b, c) # MT = M + T
        for T in Ts:
            D[-1].append(P_born(T, R, V=V) * P_adult(MT - T, I(R), γ=γ))
    return D
In [11]:
D = probability_matrix(Ts, Rs, γ=0.4, V=5.48)
fig2a_data = graphs.fig2a(D)