lectures/09_glm.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Linear Models\n",
    "\n",
    "Now that we have a better sense for neural spike traints, let's build probablistic models that predict neural responses to sensory stimuli or other covariates. These are called **encoding models**, and ideally, these models will recapitulate summary statistics of interest."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Nonlinear Poisson (LNP) models\n",
    "\n",
    "First, consider a single neuron. Let $y_{t} \\in \\mathbb{N}_0$ denote the number of spikes it fires in the $t$-th time bin. (As before, assume time bins are length $\\Delta$, typically 5-100 ms.) Let $\\mathbf{x}_t$ denote the covariates at time $t$. For example, the covariates may be features of a sensory stimulus at time bin $t$. \n",
    "\n",
    "A common modeling assumption in neuroscience is that neural spike counts are **conditionally Poisson**\n",
    "\n",
    "$$\n",
    "y_{t} \\sim \\mathrm{Po}(\\lambda(\\mathbf{x}_{1:t}) \\cdot \\Delta),\n",
    "$$\n",
    "\n",
    "where $\\mathbf{x}_{1:t} = (\\mathbf{x}_1, \\ldots, \\mathbf{x}_t)$ is the stimulus up to and including time $t$, and where $\\lambda(\\mathbf{x}_{1:t})$ is a conditional **firing rate** that depends on the stimuli.\n",
    "\n",
    "As written above, the firing rate $\\lambda$ looks like a rather complex function... it takes in an arbitrarily long stimulus history and outputs a non-negative scalar. We will make a few simplifying assumptions in order to construct our first model.\n",
    "\n",
    "1. Assume that $\\lambda$ only depends on a finite set of **features** of the stimulus history, $\\boldsymbol{\\phi}_t = (\\phi_1(\\mathbf{x}_{1:t}), \\ldots, \\phi_{D}(\\mathbf{x}_{1:t}))^\\top \\in \\mathbb{R}^D$. For example, the features may be the most recent $D$ frames of the stimulus, corresponding to $\\phi_d(\\mathbf{x}_{1:t}) = \\mathbf{x}_{t-d}$.\n",
    "\n",
    "2. Assume that $\\lambda$ only depends on **linear projections** of the features, $\\mathbf{w}^\\top \\boldsymbol{\\phi}_t \\in \\mathbb{R}$, for some weights $\\mathbf{w} \\in \\mathbb{R}^D$. We will call $\\mathbf{w}^\\top \\boldsymbol{\\phi}_t$ the **activation** at time $t$.\n",
    "\n",
    "3. Finally, assume that $\\lambda$ maps the activation through a **rectifying nonlinearity**, $f: \\mathbb{R} \\mapsto \\mathbb{R}_+$, to obtain a non-negative firing rate.\n",
    "\n",
    "Altogether, these assumptions imply a **linear nonlinear Poisson (LNP)** model,\n",
    "\n",
    "$$\n",
    "y_t \\sim \\mathrm{Po}(f(\\mathbf{w}^\\top \\boldsymbol{\\phi}_t) \\cdot \\Delta)\n",
    "$$\n",
    "\n",
    "Typical choices of rectifying nonlinearity are the exponential function, $f(a) = e^a$, and the softplus function, $f(a) = \\log (1+e^a)$. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incorporating spike history\n",
    "\n",
    "The model above treats the spike counts $y_t$ and $y_{t'}$ as **conditionally independent** given the stimulus. However, we know this assumption is invalid due to neurons' refractory period: after a neuron spikes, it cannot spike for at least a few milliseconds. For small time bins, these dependencies matter. \n",
    "\n",
    "A simple way to address this model misspecification is to allow the firing rate to depend on both the stimulus and the **spike history**, $\\lambda(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1})$. We can do so by including the spike history in the features, \n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\phi}_t = \\left(\\phi_1(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1}), \\ldots, \\phi_D(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1}) \\right)^\\top.\n",
    "$$\n",
    "\n",
    "This way, some of our features can capture the stimulus, and others can capture recent spike history. For example, one of our features might be $\\phi_d(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1}) = y_{t-d}$. In the language of statistical time series models, these spike history terms make this an **autoregressive (AR) model**.\n",
    "\n",
    ":::{admonition} Exercise\n",
    ":class: tip\n",
    "Suppose our features were $\\phi_d(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1}) = y_{t-d}$ for $d=1,\\ldots,D$. If neurons have a refractory period that prevents firing in two adjacent time bins, what would you expect the best-fitting weights $\\mathbf{w} \\in \\mathbb{R}^D$ to look like?\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-neuronal spike train models\n",
    "\n",
    "So far, we've considered models for a single neuron. In practice, we will often record from many neurons simultaneously, and we would like our models to capture correlations between neurons. \n",
    "\n",
    "Let $\\mathbf{y}_t = (y_{t,1}, \\ldots, y_{t,N})^\\top \\in \\mathbb{N}_0^N$ denote the vector of spike counts from $N$ neurons in time bin $t$. We can generalize the LNP model above as,\n",
    "\n",
    "$$\n",
    "y_{t,n} \\sim \\mathrm{Po}(f(\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t) \\cdot \\Delta)\n",
    "$$\n",
    "\n",
    "where the weights $\\mathbf{w}_n \\in \\mathbb{R}^D$ are specific to neuron $n$, and where $\\boldsymbol{\\phi}_t \\in \\mathbb{R}^D$ now includes features of the stimulus as well as the spike history of _all neurons_.\n",
    "\n",
    "For example, we might have,\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\phi}_t = (\\mathbf{x}_t,\\ldots,\\mathbf{x}_{t-L}, y_{t-1,1}, \\ldots, y_{t-L,1}, \\ldots, y_{t-1,N}, \\ldots, y_{t-L,N}, 1)\n",
    "$$\n",
    "\n",
    "where $L$ is the maximum lag of stimulus and spike history to be considered. The final 1 in $\\boldsymbol{\\phi}_t$ is a **bias** term that allows the model to learn a baseline firing rate.\n",
    "\n",
    "The entries of $\\mathbf{w}_n$ associated with the features $(y_{t-1,m}, \\ldots, y_{t-L,m})^\\top$ can be thought of as **coupling filters**, which model how spikes on neuron $m$ influence the future firing rate of neuron $n$. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis function encodings\n",
    "\n",
    "The model above has $\\mathcal{O}(N^2 L)$ weights for the coupling filters. For small bin sizes, $L$ may need to include dozens of past time bins to capture all the pairwise interactions. However, these coupling filters are often approximately smooth functions of the time lag. One way to cut down on parameters and capture this smoothness is to use a **basis function representation**. For example, one of the features can be,\n",
    "\n",
    "$$\n",
    "\\phi_{m,b}(\\mathbf{x}_{1:t}, \\mathbf{y}_{1:t-1}) = \\sum_{\\ell=1}^L y_{t-\\ell,m} e^{-\\frac{1}{2 \\sigma^2}(\\ell - \\mu_b)^2}.\n",
    "$$\n",
    "\n",
    "This is a **radial basis function** encoding of the spike history of neuron $m$. It is a weighted sum of past spiking, where the weights are a squared exponential (aka Gaussian) kernel centered on delay $\\mu_b$. We can use $B < L$ basis functions to summarize the spike history over the last $L$ time bins.\n",
    "\n",
    ":::{admonition} Exercise\n",
    ":class: tip\n",
    "Show that the feature above can be written as a convolution of the spike history with a squared exponential kernel.\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalized linear models (GLMs)\n",
    "\n",
    "The model described above, with stimulus features, spike history terms, and basis function encodings, is what neuroscientists often call \"the\" generalized linear model (GLM), after {cite:t}`pillow2008spatio`. Of course, in statistics we know that this model is just one instance of a broad family of GLMs, which are characterized by linear projections of covariates, nonlinear link functions, and exponential family conditional distributions {cite}`mccullagh2019generalized`. In fact, we have already encountered one GLM in this course: the logistic regression model from [Unit 1](./07_pose_tracking.ipynb)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum likelihood estimation (MLE)\n",
    "\n",
    "As with logistic regression, the GLMs constructed above won't have closed form solutions for the MLE. However, they will be amenable to the optimization methods discussed previously: (stochastic) gradient descent, Newton's method, etc.\n",
    "\n",
    "Let's follow the same recipe as before to compute the gradient of the negative log likelihood of a GLM with an exponential nonlinearity,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{L}(\\mathbf{w}) \n",
    "&= - \\log p(\\mathbf{y} \\mid \\mathbf{x}, \\mathbf{w}) \\\\\n",
    "&= -\\sum_{t=1}^T \\sum_{n=1}^N \\log \\mathrm{Po}(y_{t,n} \\mid e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta) \\\\\n",
    "&= -\\sum_{t=1}^T \\sum_{n=1}^N \\left( -\\log y_{t,n}! + y_{t,n} \\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t  + y_{t,n} \\log \\Delta - e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta \\right)  \\\\\n",
    "&= -\\sum_{t=1}^T \\sum_{n=1}^N \\left(y_{t,n} \\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t - e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta \\right)  + c\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $c$ is constant wrt $\\mathbf{w}_n$.\n",
    "\n",
    "```{admonition} Separability\n",
    "The loss function is a sum over neurons, and that the weights for neuron $n$ do not interact with those of other neurons $m$. That means we can write the objective function as,\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\mathbf{w}) = \\sum_{n=1}^N \\mathcal{L}_n(\\mathbf{w}_n),\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\mathcal{L}_n(\\mathbf{w}_n) = - \\left(y_{t,n} \\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t - e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta \\right).\n",
    "$$\n",
    "\n",
    "We can optimize the weights for each neuron separately. (And in parallel, if you like!)\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradients and Hessians\n",
    "\n",
    "Now take the gradient wrt the weights of neuron $n$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla_{\\mathbf{w}_n} \\mathcal{L}_n(\\mathbf{w}_n) \n",
    "&= -\\sum_{t=1}^T (y_{t,n} - e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta) \\boldsymbol{\\phi}_t \\\\\n",
    "&= -\\sum_{t=1}^T (y_{t,n} - \\mathbb{E}[y_{t,n} \\mid \\mathbf{w}_n] ) \\boldsymbol{\\phi}_t \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The gradient takes an intuitive form: it is a sum of the feature vectors, $\\boldsymbol{\\phi}_t$, weighted by the **error** $(y_{t,n} - \\mathbb{E}[y_{t,n} \\mid \\mathbf{w}_n])$. Recall that we found the same simple form  for the gradient of the logistic regression model. We'll show that this form is characteristic of a special family of GLMs in just a second.\n",
    "\n",
    "First, note that the negative log likelihood is again a convex function. We can see that from the Hessian,\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}_n}^2 \\mathcal{L}_n(\\mathbf{w}_n) \n",
    "= \\sum_{t=1}^T (e^{\\mathbf{w}_n^\\top \\boldsymbol{\\phi}_t} \\cdot \\Delta) \\, \\boldsymbol{\\phi}_t \\boldsymbol{\\phi}_t^\\top.\n",
    "$$\n",
    "\n",
    "Since this is a weighted sum of outer products with positive weights, the Hessian is positive semi-definite (PSD) and the negative log likelihood is convex. Convexity implies that gradient descent and other optimization methods will (under weak conditions) find a global optimum.\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential family distributions\n",
    "\n",
    "The reason why the Poisson GLM and logistic regression had such nice mathematical forms is because they are both **exponential family** GLMs. Exponential family distributions have densities of the form,\n",
    "\n",
    "$$\n",
    "p(y; \\theta) = h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle - A(\\theta) \\right\\}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "- $\\theta$ is the **natural parameter** of the distribution\n",
    "- $t(y)$ is a **sufficient statistic** of the data\n",
    "- $h(y)$ is the **base measure**, which we won't worry too much about\n",
    "- $A(\\theta)$ is the **log normalizer**, which ensures the density integrates to one.\n",
    "\n",
    "In order for the density to integrate to one, the log normalizer must be,\n",
    "\n",
    "$$\n",
    "A(\\theta) = \\log \\int h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle \\right\\} \\, \\mathrm{d} y\n",
    "$$\n",
    "\n",
    "Many of the distributions we've encountered thus far are exponential family distributions: Poisson, Bernoulli, Gaussian, gamma, etc. \n",
    "\n",
    "```{admonition} Example\n",
    "The Poisson distribution can be written in exponential family form by rearranging its pmf to,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathrm{Po}(y; \\lambda) \n",
    "&= \\frac{1}{y!} \\lambda^y e^{-\\lambda} \\\\\n",
    "&= \\frac{1}{y!} \\exp \\left\\{ y \\log \\lambda - \\lambda \\right\\} \\\\\n",
    "&= \\frac{1}{y!} \\exp \\left\\{ y \\theta - e^\\theta \\right\\} \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $t(y) = y$ is the sufficient statistic, $\\theta = \\log \\lambda$ is the natural parameter, $e^\\theta$ is the log normalizer, and $h(y) = 1/y!$ is the base measure.\n",
    "```\n",
    "\n",
    "```{admonition} Exercise\n",
    ":class: tip\n",
    "Write the Bernoulli pmf, \n",
    "\n",
    "$$\n",
    "\\mathrm{Bern}(y; p) = p^y (1-p)^{1-y},\n",
    "$$\n",
    "\n",
    "in exponential family form. What is the natural parameter, the sufficient statistic, the base measure, and the log normalizer?\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential family GLMs \n",
    "\n",
    "An exponential family GLM is one in which the observation distribution is an exponential family and the natural parameter is a linear-nonlinear function of the features,\n",
    "\n",
    "$$\n",
    "\\theta = g(\\mathbf{w}^\\top \\boldsymbol{\\phi}),\n",
    "$$\n",
    "\n",
    "where $g: \\mathbb{R} \\mapsto \\mathbb{R}$ is a monotonically increasing function.\n",
    "\n",
    "<!-- The canonical choice is to make $g$ the identity function.\n",
    "\n",
    "Buf if we make $g$ the identiy function, then we simply have a linear model. Where did the _generalized_ (i.e., nonlinear) part go? Note that the model is linear in the **natural parameters** $\\theta$ of the distribution, whereas the Poisson GLM and logistic regression were defined in terms of the **mean parameters**. -->\n",
    "\n",
    "Note that this model is defined in terms of the natural parameters $\\theta$ rather than the mean parameters $\\mathbb{E}[t(y) \\mid \\theta]$ (which for a Poisson distribution is simply $\\mathbb{E}[y \\mid \\theta]$).  How do the natural and mean parameters relate to one another? Through the gradient of the log normalizer,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "A'(\\theta) &\\triangleq \\frac{\\mathrm{d}}{\\mathrm{d} \\theta} A(\\theta) \\\\\n",
    "&= \\frac{\\mathrm{d}}{\\mathrm{d} \\theta} \\log \\int h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle\\right\\} \\, \\mathrm{d}y \\\\\n",
    "&= \\frac{\\int t(y)  h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle\\right\\} \\, \\mathrm{d}y}{\\int h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle\\right\\} \\, \\mathrm{d}y} \\\\\n",
    "&= \\int t(y)  h(y) \\exp \\left\\{ \\langle t(y), \\theta \\rangle - A(\\theta) \\right\\} \\, \\mathrm{d}y \\\\\n",
    "&= \\mathbb{E}[t(y) \\mid \\theta]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "In short, _gradients of the log normalizer yield expected sufficient statistics!_\n",
    "\n",
    "<!-- When the sufficient statistics are linearly independent, the exponential family is said to be **minimal**, and there is a 1:1 mapping between the natural and mean parameters. That is, in a minimal exponential family, the gradient of the log normalizer is invertible. -->\n",
    "\n",
    "Now let's write the GLM in terms of the mapping to the mean parameters,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbb{E}[t(y) \\mid \\mathbf{w}, \\boldsymbol{\\phi}]\n",
    "= A'(\\theta) \n",
    "= f(\\mathbf{w}^\\top \\boldsymbol{\\phi})\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $f(a) = A'(g(a))$ is the **mean function**. The **canonical mean function** is obtained by setting $g$ to the identity function, so that $f$ is simply the gradient of the log normalizer. \n",
    "\n",
    "For example, in the Poisson distribution, the sufficient statistic is $t(y) = y$, and the log normalizer is $A(\\theta) = e^\\theta$. The derivative of the log normalizer is also $A'(\\theta) = e^\\theta$, so the canonical Poisson GLM is,\n",
    "\n",
    "$$\n",
    "\\mathbb{E}[y \\mid \\mathbf{w}, \\boldsymbol{\\phi}] = e^{\\mathbf{w}^\\top \\boldsymbol{\\phi}_t},\n",
    "$$\n",
    "\n",
    "just like we had above! (Note: the equivalence is exact when bin size is $\\Delta = 1$; the scaling factor doesn't change the story much.)\n",
    "\n",
    "<!-- For example, in the Poisson GLM described above, we modeled the Poisson mean as $\\lambda = e^{\\mathbf{w}^\\top \\boldsymbol{\\phi}}$. (Assume a bin size of $\\Delta = 1$ for simplicity.) To put this in exponential family form, we need the natural parameter  -->\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradients for canonical exponential family GLMs\n",
    "\n",
    "Now consider the negative log likelihood in a canonical exponential family GLM,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{L}(\\mathbf{w})\n",
    "&= - \\sum_{i=1}^T \\log p(y_i \\mid \\mathbf{w}, \\boldsymbol{\\phi}_i) \\\\\n",
    "&= - \\sum_{i=1}^T \\langle t(y_i), \\mathbf{w}^\\top \\boldsymbol{\\phi}_i \\rangle - A(\\mathbf{w}^\\top \\boldsymbol{\\phi}_i) + c\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The gradient is,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla \\mathcal{L}(\\mathbf{w})\n",
    "&= - \\sum_{i=1}^T (t(y_i) - A'(\\mathbf{w}^\\top \\boldsymbol{\\phi}_i)) \\, \\boldsymbol{\\phi}_i \\\\\n",
    "&= - \\sum_{i=1}^T (t(y_i) - \\mathbb{E}[t(y_i) \\mid \\mathbf{w}, \\boldsymbol{\\phi}_i]) \\, \\boldsymbol{\\phi}_i.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Thus, we see that the simple form of the gradient &mdash; a sum of features weighted by the corresponding error &mdash; follows from a general property of exponential family GLMs.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLMs in neuroscience\n",
    "\n",
    "Generalized linear models are central tools in systems neuroscience. Work by {cite:t}`paninski2004maximum`, {cite:t}`truccolo2005point`, and {cite:t}`pillow2008spatio` laid the theoretical foundations for much of the experimental work that followed. \n",
    "\n",
    "- {cite:t}`park2014encoding`, {cite:t}`yates2017functional` applied GLMs to characterize heterogeneous responses in multi-neuronal spike train recordings from lateral intraparietal cortex (LIP) during decision making. \n",
    "\n",
    "- {cite:t}`hardcastle2017multiplexed` applied a similar analysis to characterize functional cell types in medial entorhinal cortex (MEC), finding grid cells, border cells, and much more.\n",
    "\n",
    "- {cite:t}`vidne2012modeling` characterized the responses of retinal ganglion cells, allowing for correlated noise processes, and {cite:t}`freeman2015mapping` further refined the model to include common inputs from bipolar cells. {cite:t}`mcintosh2016deep` framed these models as convolutional neural networks, which we will study in lab.\n",
    "\n",
    "- On the methodological side, much work has gone into fitting these models at scale {cite}`ramirez2014fast,zoltowski2018scaling`, extending the models with structured prior distributions {cite}`linderman2016bayesian`, making them more biophysically realistic {cite}`latimer2019inferring`, and understanding the extent of spiking responses that GLMs can capture {cite}`weber2017capturing`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.7.11 (default, Jul 27 2021, 07:03:16) \n[Clang 10.0.0 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5b28c5bd4ee93d765ebe901023d5522822fb8ad083dac3187c5545022f913719"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}