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 routines, using the bokeh library, have been abstracted away in the graphs.py file.
A citable version of this notebook is available at figshare. You can contact me for questions or remarks at
import numpy as np import graphs np.random.seed(0) # repeatable results.
Piantosi and Kidd's model ties together three main parameters:
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.
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) when the head size at birth exceed a fixed parameter $V$, fitted to 5.48 cm. Here $T$ is the duration of the gestation period.
def phi(z): """Sigmoid function""" return 1/(1 + np.exp(-z)) def P_sb(T, R, V=5.48, g=g): """Probability to survive birth""" return phi(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$.
Here, the parameter $\gamma$ is fitted to 0.4.
def P_sc(M, I_p, gamma=0.4): """Probability to survive childhood""" return np.exp(-max(0, M) * gamma / I_p)
The article assumes that $I_p$ is equal to $R$  for Figure 1 and 2A. To try other relationships, modify the following function:
def I(R): """Return the intelligence (of the parents) as a function of R""" return R
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.
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_sb, P_sc, I, R=8.4/2)
def update_fig1(R): graphs.update_fig1(fig1_data, g, P_sb, P_sc, 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.
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$. We employ a simple dichotomy method regardless, rather than a precomputed answer, to allow for arbitrary modification of the $g$ function.
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. It may not hold if you modify the g function # above. In this case, comment this line, the code below is general. 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)
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, gamma=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_sb(T, R, V=V) * P_sc(MT - T, I(R), gamma=gamma)) return D
D = probability_matrix(Ts, Rs, gamma=0.4, V=5.48) fig2a_data = graphs.fig2a(D)