Skip to content

Code for training and validation of ion channel models fitted to sine-wave based voltage clamp data

License

Notifications You must be signed in to change notification settings

mirams/sine-wave

Repository files navigation

Sinusoidal voltage protocols for rapid characterization of ion channel kinetics: Supplementary Data

This is code associated with the paper:
Beattie, K. A., Hill, A. P., Bardenet, R., Cui, Y., Vandenberg, J. I., Gavaghan, D. J., de Boer, T. P. & Mirams, G. R. Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics. J. Physiol. 596, 1813–1828 (2018).

In this readme we detail the contents of this Supplementary Data repository. We'll explain how all the scripts in this folder are used to generate the results in the paper, and which scripts are used to collate and plot the data for each figure.

Table of Contents

In Protocols we include the voltage clamp waveform for each of the protocols. These files include a list of voltages that comprise the protocol (in mV), with each row corresponding to a 0.1 ms timestep (10kHz samples).

The processed (leak and dofefilide subtracted) experimental data are included in ExperimentalData. There is a folder corresponding to the data for each cell (the correspondence between file names and cell numbers in the paper is provided in cell_index.txt). There are additional folders containing the average data for the sine wave ('average') and one containing the repeats from the sine wave protocol for each cell ('sine_wave'). These data traces correspond to leak and dofetilide subtracted data.

Links to the full set of raw data traces in both .abf and plain text format can be found in FullExperimentalData.

Here we provide details of all the codes for

  1. Running simulations with a given model and parameter set
  2. Minimisation (CMA-ES)
  3. MCMC
  4. Plotting the figures presented in the manuscript

You should set up the below programs before attempting to run our code.

The parameter values for each model are included in ParameterSets. Note that in each parameter set the final parameter is the conductance parameter which has been set to 0.1 for all models. This value is irrelevant as we scale the literature model simulations to either a simulated or experimental reference trace when plotting these model simulations or using them for comparison.

For each model (or sets of models which share the same model structure and only vary in their parameterisations) there is a Mex file detailing the set of equations which define that model. For example MexWang.c is used to simulate the Wang et al. 1997 model. The appropriate Mex file is called when SimulatingData.m is run.

If any changes to Mex files are made these must be recomplied using:

mex -I/path_to_CVODE/include -L/path_to_CVODE/lib -lsundials_nvecserial -lsundials_cvode MexFilename.c

modeldata.m defines the model_type for each model (to determine which Mex file should be used for each model simulation) and also identifies the appropriate parameter set in ParameterSets to be used when simulating each model.

To find the best fit to the sine wave experimental data we first run FullGlobalSearch.m for each cell and then once we've verified that the CMA-ES algorithm repeatedly returns parameters in the same region of parameter space on multiple different iterations we then run AdaptiveMCMCStartingBestCMAES.m to determine MCMC chains.

cmaes.m defines the CMA-ES algorithm used for the initial search of the parameter space before running MCMC. This file was downloaded from https://www.lri.fr/~hansen/cmaes_inmatlab.html

FullGlobalSearch.m is run to search the parameter space and identify the best guess parameter set to be used to initialise the MCMC chain.

AdaptiveMCMCStartingBestCMAES.m defines the covariance adaptive MCMC algorithm used.

The contents of CMAESFullSearchResults give an example output of CMAES search results to be used for running of the MCMC algorithm (for cell #5).

In MCMCResults there are all the MCMC chains run from fitting the model to each cell's sine wave data (as well as the average model). We include a separate file for each cell for: (i) the parameter values in the MCMC chain; (ii) the likelihoods corresponding to these parameter values; and (iii) the acceptance rate over the running of the MCMC chain.

PlottingSamplesFromMCMC.m is used to assess the 95% credible intervals for the fits to the sine wave protocol when taking 1000 samples from the MCMC parameter distribution results. This figure is quoted in the text in Model Calibration Section 2.2.

In SimulatedData we have the simulated data trace from the maximum likelihood parameters identified from fitting to experimental data for cell 5. This data trace was produced using ProducingSimulatedDataWithNoise.m. We use these data to produce the simulated data trace MCMC distributions shown in Figure C5.

In MCMCResultsSimulated there is the MCMC chain (and likelihood and acceptancerate files) for the synthetic data study results shown in Figure C5.

CreatingAveragedModel.m creates the 'averaged' sine wave data which was used to fit an 'average' model for comparison of cell-specific vs. average model predictions.

Code to generate and plot data for each of the figures in the main text and supplement are listed here:

Any remaining .m files not described in this document are used within the individual scripts described above with a note within each script about their use.

Note that the Zeng model contained singularities and was unable to simulate our action potential protocol (Pr6) in its original form, so whenever the action potential waveform is simulated we use a version of the Zeng et al. model where the singularities have been removed (as in MexZengTol.c).

We also note that the leak subtracted dofetilide subtracted experimental traces for the activation kinetics protocols for cells 3, 4 and 7 (16708016, 16708060 and 16704007) appear to be over leak subtracted when we use the same leak resistance values identified for appropriately leak subtracting the sine wave trace for these cells. The activation kinetics protocols are performed at the very start of the protocol sequence so we may expect some change in leak resitance by the time the sine wave protocol is performed. As we do not use the activation kinetics data for these cells and the remaining protocols appear to be appropriately leak subtracted using the leak resistance values identified for each cell from the sine wave protocol, we do not adjust the leak resistance for activation kinetics data.

A final point to note is that when comparing our input voltage protocol with the command voltage recorded in the results file output from the pClamp10 software, we observed a time shift of one sample period (0.1 ms) in the applied protocol to that recorded in experiment. This shift was consistent in all recordings, and is likely due to a timing offset in the digitizer. To account for this observation, to ensure the use of recorded currents and correct corresponding applied voltage, we applied a time shift of 0.1 ms to the input protocol when performing simulations for comparison with experimental data.