Analysis of activity dependent development of brain topographic maps​
​by Nicholas Gale

Article Abstract

Topographic maps are a brain structure connecting pre-synpatic and post-synaptic brain regions. Topographic development is dependent on Hebbian-based plasticity mechanisms working in conjunction with spontaneous patterns of neural activity generated in the pre-synaptic regions. Studies performed in mouse have shown that these spontaneous patterns can exhibit complex spatial-temporal structures which existing models cannot incorporate. Neural field theories are appropriate modelling paradigms for topographic systems due to the dense nature of the connections between regions and can be augmented with a plasticity rule general enough to capture complex time-varying structures. We propose a theoretical framework for studying the development of topography in the context of complex spatial-temporal activity fed-forward from the pre-synaptic to post-synaptic regions. Analysis of the model leads to an analytic solution corroborating the conclusion that activity can drive the refinement of topographic projections. The analysis also suggests that biological noise is used in the development of topography to stabilise the dynamics. MCMC simulations are used to analyse and understand the differences in topographic refinement between wild-type and the β2 knock-out mutant in mice. The time scale of the synaptic plasticity window is estimated as 0.56 seconds in this context with a model fit of R2=0.81.

Overview

This notebook acts as a companion to and generates the analysis of the paper: “Analysis of Activity Dependent Development of Topographic Maps in Neural Field Theory with Short Time Scale Dependent Plasticity”
1
​
. In brief, we are interested in the role of Hebbian plasticity and spontaneously generated activity in the formation of a ubiquitous biological system: the topographic map. Current modelling of the role activity plays in the development of such systems is limited to assumptions of static activity input; a region of the receptive field of the post-synaptic region will be stimulated until equilibrium. The data generated by imaging experiments in mice show that the spontaneously generated activity comes in the form of propagating waves with complex spatio-temporal statistics. We then develop a framework in which the transient nature of these waves can be incorporated. In this notebook we will analyse the consequences and predictions of the model and examine how well the model can explain a key organisational feature of topographic maps: the width of the synaptic kernel. This is done by Markov-Chain-Monte-Carlo analysis.

Packages

This analysis will rely on a few dependencies coming both from Mathematica packages and a piece of custom software written by Josh Burkart to perform MCMC sampling
2
​
. A slight modification to the package has been made to deal with boundary conditions that are not at infinity because we will be using uniform priors that terminate the sample space. The modification prevents the MCMC from attempting to sample outside of this space. We will load the native packages into Mathematica’s parallel framework but the MCMC functions will be imported to the parallel threads at a later time when those threads are initialized.
Needs["FourierSeries`"]​​Needs["HypothesisTesting`"]​​ParallelNeeds["FourierSeries`"]​​ParallelNeeds["HypothesisTesting`"]​​Get[FileNameJoin[{NotebookDirectory[],"mcmc.m"}]];​​$Messages={};

Model

Here we will define the deterministic model derived from the analysis presented in the paper. This notebook will not detail the precise derivations that lead to the model; it is mainly concerned with the computational aspects of the paper and the generation of statistics which encapsulate biological descriptions and predictions. In brief, the model details the development of topographic maps which, in one dimension, can be characterised by a fundamental length referred to here as width. The model is developed quite generally but, as with all models, is not scientifically useful without comparison to data. The data we aim to predict comes from wild-type (genetically unaltered) and β2 (genetically altered) mice. The genetic alteration disrupts neural activity which the model is able to capture, and the phenotype of genetically altered mice increases the characteristic width of the organisation making these two references useful benchmarks for the model. By calibrating against the measurements we are able to make predictions about the effect of various biological mechanisms by interrogating the model.
​
Here we will also instantiate the common or default parameters that won' t change the estimate for the width of the topographic organisation which are primarily interested in. We should also keep in mind that while these parameters will not affect the prediction for the width they certainly affect other features of the kernel - the complete description of synaptic connections in the system - and we do not claim that these are necessarily the correct parameters for our biological system of interest, simply that they will not affect the analysis.

Deterministic Model

(*Firingrate*)​​Q[rate_,v_,β_,t_]:=rate/(1+Exp[-β*(v-t)])​​​​(*LateralConnectivity*)​​W[x_,r1_,r2_,R1_]:=R1*r1Exp[-1/2*x^2r1^2]-r2Exp[-1/2*x^2r2^2];​​​​(*Wave-Form*)​​h[k_,s1_,s2_]:=s1s2*Exp[-1/2*k^2*(s2)^2]​​​​(*Non-ConstantPartofFT{S}(k)*)​​g[k_,c_,τ_,tp_,r1_,r2_,R2_]:=4c^2k^2τAbs[tp]/(1+c^2k^2tp^2)/(c^2τ^2k^2+(1-W[k,r1,r2,R2])^2)​​​​(*Fouriertransformofthesynaptickernel*)​​ftS[k_,κ_,λ_,τ_,τprime_,tp_,ap_,fr_,β_,θ_,s1_,s2_,c_,r1_,r2_,R1_,u0_]:=Re[κ/λ(1/(1-g[k,c,τ,tp,r1,r2,R1]*h[k,s1,s2]^2*ap^2*Exp[θ*β-β*u0]*fr*β/(1+Exp[θ*β-β*u0])^2*τprime/λ)-1)]​​ftSrenorm[k_,κ_,λ_,τ_,τprime_,tp_,ap_,fr_,β_,θ_,s1_,s2_,c_,r1_,r2_,R1_,u0_]:=Re[κ/λ(1/(1-g[k,c,τ,tp,r1,r2,R1]*h[k,s1,s2]^2*ap^2*Exp[θ*β-β*u0]*fr*β/(1+Exp[θ*β-β*u0])^2*τprime/λ)-1)]-κ/λ​​​​(*SynapticKernel*)​​S[x_,κ_,λ_,τ_,τprime_,tp_,ap_,fr_,β_,θ_,s1_,s2_,c_,r1_,r2_,R1_,u0_]:=Re[NInverseFourierTransform[ftS[k,κ,λ,τ,τprime,tp,ap,fr,β,θ,s1,s2,c,r1,r2,R1,u0],k,x,Method"GaussKronrodRule"]]

Parameters

Here we define default parameters which are used throughout the paper, and this notebook, unless otherwise specified.
κ=0.001;​​λ=0.001;​​ap=1;​​β=0.26;​​θ=-45;​​s1=0.05;​​fr=1;​​τprime=100;​​u0=-58;

Analysis Functions

Here we define two functions which will be useful for analysis: the width of the organisation, and the stability function. The organisation width is the predictor for the biological measurements made in the mouse superior colliculus retinotopic map. The stability function determines the points at which the model begins to exhibit singularities and therefore becomes inappropriate to explain biological data.
(*Widthestimate*)​​SWidth[τ_?NumericQ,tp_?NumericQ,s2_?NumericQ,c_?NumericQ,r1_?NumericQ,r2_?NumericQ,R1_?NumericQ]:=Abs[1/(10^-7+ArgMax[ftS[k,κ,λ,τ,τprime,tp,ap,fr,β,θ,s1,s2,c,r1,r2,R1,u0],k])]​​(*SWidth[τ_?NumericQ,tp_?#<=0&,s2_?NumericQ,c_?NumericQ,r1_?NumericQ,r2_?NumericQ,R1_?NumericQ]:=10^9*)​​​​(*Stabilityfunction*)​​stab[k_,τ_?NumericQ,tp_?NumericQ,s2_?NumericQ,c_?NumericQ,r1_?NumericQ,r2_?NumericQ,R1_?NumericQ]:=g[k,c,τ,tp,r1,r2,R1]*h[k,s1,s2]^2*ap^2*Exp[θ*β-β*u0]*fr*β/(1+Exp[θ*β-β*u0])^2*τprime/λ​​stabMax[τ_?NumericQ,tp_?NumericQ,s2_?NumericQ,c_?NumericQ,r1_?NumericQ,r2_?NumericQ,R1_?NumericQ]:=Max[stab[Table[k,{k,0.0,10,0.01}],τ,tp,s2,c,r1,r2,R1]]

Parameter Sensitivity

The model has several parameters which have various effects on model behaviour and predictions. A key predictor is the organisation width defined in the previous section which is sensitive only to a few parameters: it is invariant to all parameters chosen expect those in
∂g
∂k
. The following widget allows informal exploration of the these parameters dependencies.
Manipulate[Plot[ftS[k,κ,λ,τ,τprime,tp,ap,fr,β,θ,s1,s2,c,r1,r2,R1,u0],{k,-1,1}],{{κ,1},0,10},{{λ,1},0,10},{{τ,1},0,10},{{τprime,1},0,10},{{tp,1},0,10},{{ap,1},0,10},{{fr,1},0,10},{{β,1},0,10},{{θ,1},-50,10},{{s1,1},0,10},{{s2,1},0,10},{{c,1},0,10},{{r1,1},0,10},{{r2,1},0,10},{{R1,1},0,10},{{u0,0},-55,0},SaveDefinitionsTrue]
​
κ
λ
τ
τprime
tp
ap
fr
β
θ
s1
s2
c
r1
r2
R1
u0
-1.0
-0.5
0.5
1.0
0.02
0.04
0.06
0.08
0.10
0.12

Statistical Analysis

We would like to fit our model to a sample of wild-type mice and β2 mice. A reliable measurement of wave properties comes from Stafford et al (2009)
3
​
. It is most appropriate to choose a measurement error model with categories of β2 and Wild-Type. In this section we are going to use an MCMC fitted model to regress on the data. We are choosing seconds as the time scale, and mm as the length scale. We will also need to estimate the parameters which dictate the recurrent kernel as this is not well studied.

Estimation of Recurrent Kernel

An essential feature of the model is the lateral or recurrent connections W. In mice there have been measurements of the extent to which the inhibitory and excitatory populations extend lengthwise in the superior colliculus; the post-synaptic target of interesting. We can estimate the data from the results shown in Figure 1 from Phongphanphanee et. al. (2014)
4
​
. These data come in the form {location (μm), Maximal PSP}. We then perform a non-linear model fit on the ansatz form of these lateral connections commonly taken in the literature: the difference of two Gaussians. This informs a “best initial guess” of the parameters but we shall also estimate them via MCMC.
Fig1LM={{-450,-17.5},{-300,-10},{-150,-20},{0,80},{150,-5},{300,-25},{450,0}};​​Fig1RC={{-450,-1},{-300,1},{-150,40},{0,85},{150,5},{300,-30},{450,10}};​​dataAv=(Fig1LM+Fig1RC)/2;​​data=Table[{dataAv[[i]][[1]]/1000,dataAv[[i]][[2]]/1000},{i,1,Length[dataAv]}];​​model[x_,R1_,r1_,r2_]:=R1*Exp[-x^2/(2r1^2)]-Exp[-x^2/(2r2^2)]​​fit=NonlinearModelFit[data,model[x,R1est,r1est,r2est],{R1est,r1est,r2est},x]​​fit["ParameterConfidenceIntervalTable"]​​Plot[fit[x],{x,-0.5,0.5}]
FittedModel
1.08248
-30.082
2
x

-
-26.9677
2
x


Estimate
Standard Error
Confidence Interval
R1est
1.08248
0.010192
{1.05418,1.11078}
r1est
0.128923
0.012295
{0.094787,0.16306}
r2est
0.136164
0.0136187
{0.0983525,0.173976}
-0.4
-0.2
0.2
0.4
-0.02
0.02
0.04
0.06
0.08

Data

We now do the necessary unit conversions to be appropriate for the model. We assume that the wave speed is measured in mm and therefore the area measurements should be scaled by a factor of ten so that s2 is also on the scale of millimeters.
We can now define the likelihood function given a data point is from a particular category and under the assumption that the other parameters are independent of these particular mutations. Most of the parameters do not affect the size and can therefore be excluded, some of the parameters can be estimated (r1, r2) while others we attach an uninformative prior (τ, tp, R1). We start by defining the deterministic model, priors, and categorical likelihood. First make a convenience wrapper for all of the data.

Priors

Need to restrict these to the positive domain which is a problem for the MCMC package. Create a new PDF and use the HeavisideTheta function to get a functional evaluatable form for the MCMC package and then add a small perturbation to avoid a simplification by the WolframKernel.

Likelihood Model

There are quite a few parameters so we shall assemble the likelihood in the following form: independent variables sampled from measurement distributions, likelihood of dependent variable under normal error assumption, and priors for parameters.
Now we define the log-likelihood for all of the data. This is what we will be minimising via MCMC.

MCMC

Now we perform the MCMC giving it lots of room to breath. We need to check that the MCMC did converge which probably will entail running several chains and doing some convergence measures. We will do 8 chains corresponding to two initializations in each of the sample basins and run for 100000 trials. Once each of the chains are complete we compute the Rubin-Gelman statistic to see if the MCMC has converged and combine all of the parameter runs together. Then we can plot the distributions, time-series, histogram and pick most likely parameters on all parameters that haven’t been estimated from data. Whether any of these other parameters in the model correlate with these experimental categories remains to be seen, but as a reasonable judgement we are assuming that they are independent: β2 knock-out is probably not going to change the time scale on which plasticity occurs, for example.

Generate the Gelman-Rubin convergence diagnostic

Goodness of Fit

We need to establish a goodness of fit for our model. While there are always a number of potential issues concerning model fits for our purposes we shall simply calculate the R2 coefficient.

Image Generation

In this section we want to generate the images that guide the paper. We will define some global plotting parameters and in each subsection highlight any plot specific parameters.

Model Summary

This is the hardest figure to do programmatically as it required many changes on the basis of input from several different people. As a result the most straightforward, albeit not very reproducible, method to create the image was to generate four individual graphics and use the Mathematica image editing tool to arrange them. The graphics consist of a base summary image of the model and three additional images containing information about various aspects of the model which are overlain over the base image. The code for generating the separate graphics is shown below and the final edited image is presented in a single cell.

Retinotopy

Typical Organisations

Width, and Amplitudes

Posteriors

References

1. Gale, N., Rodger, J., Small, M. & Eglen, S. Analysis of activity dependent development of topographic maps in neural field theory with short time scale dependent plasticity. Mathematical Neuroscience and Applications 2, (2022)
2. J. Burkart. Mathematica Markov chain Monte Carlo. https://github.com/joshburkart/mathematica-mcmc, (2017)
3. Stafford, B. K., Sher, A., Litke, A. M. & Feldheim, D. A. Spatial-temporal patterns of retinal waves underlying activity-dependent refinement of retinofugal projections. Neuron 64, 200–212 (2009)
4. Phongphanphanee, P. et al. Distinct local circuit properties of the superficial and intermediate layers of the rodent superior colliculus. Eur. J. Neurosci. 40, 2329–2343 (2014)