Appendix

The appendix includes additional simulation experiments, detailed implementation of the algorithm, stimulus reconstruction analysis, goodness-of-fit tests, and a discussion about the inverse mapping.

1.1 A. Simulation study

First, we created a series of PP-GLMs with coefficients \(\varvec^0(g), g \in \\). These parameters corresponded to the models with different ion channel conductance scaling factors. Differences between adjacent models \(\varvec^0(g_i)\) and \(\varvec^0(g_)\) were small and the trend was smooth. The parameters came from a previous fit. Then for each \(\varvec^0(g_i)\), we simulated 100 3-second spike trains according to Eq. (1). The simulated spike trains were then used to fit new PP-GLMs in Eq. (6), and the penalty hyperparameter \(\lambda\) was selected as described above in Eq. (7). We expected to see that after applying the trend filtering technique, the model could recover the trend of changes despite the Poisson-like noise from the spike trains. We repeated the above procedures 100 times to acquire the mean and the variance of the error. Besides trend recovery simulation, we also checked the goodness-of-fit using KS test based on time rescaling theorem (Brown et al., 2002; Haslinger et al., 2010). All fitted models had good performance (data not shown). The results are shown in Fig. 6.

Fig. 6

Simulations verification for the joint training model 6. A, B, and C provide one example fit. D summarizes 100 repeated fits. SS values for PP-GLM fits of simulated spike trains for (A) the stimulus and (B) post-spike history coefficients, and C the log-likelihood all as a function of \(\lambda\). True SS values are shown at the bottom of A and B. Panels are similar to Fig. 3K-M, except that the true model refers to a known set of PP-GLMs with coefficients \(\varvec^0(g), g \in \\). When \(\lambda = \lambda ^*\) (gray dashed line), the SS values are very close to the true SS values, thereby validating our trend filtering penalty hyperparameter selection method. D SS error between true PP-GLM and 100 separate sets of simulated spike trains from the true PP-GLM as a function of aligned \(\lambda\) index. Since different runs may choose different optimal tuning parameter, so the tuning parameters along the x-axis, the index \(\lambda\), are aligned to the optimal \(\lambda ^*\) at index 0. The index in the plot is normalized by \(\lambda - \lambda ^*\)

Fig. 7

Sensitivity analysis. The histogram presents the distribution of selected hyperparameter \(\lambda _*\) with different threshold \(\zeta\) (as in Eq. (7)) using the same set of simulations as in Fig. 6. The penalty \(\lambda\) is selected among a set of discrete values \(\lambda \in \Lambda = \, \lambda _ \alpha , \lambda _ \alpha ^2,... ,\lambda _\alpha ^, 0 \}\), where \(k = 22\) and \(\alpha = e^\)

Figure 7 shows the distribution of selected hyperparameter \(\lambda _*\) with different threshold \(\zeta\) using the same set of simulations as in Fig. 6. The selection method is in Eq. (7). When \(\log 1.0001 \le \zeta \le \log 1.001\), the distribution of \(\lambda _*\) almost remains the same. When \(\zeta \ge \log 1.05\) is a large value, the selected \(\lambda _*\) is shifted toward right a little, but the value of \(\lambda _*\) does not become much larger. This is because when the penalty becomes stronger, the performance or the log-likelihood drops dramatically. This implies the hyperparameter selection method is not sensitive to the threshold \(\zeta\) between \(\log 1.0001\) and \(\log 1.05\).

1.2 B. ADMM optimization algorithm for training PP-GLMs with trend filtering1.2.1 B.1. Update rules

Training PP-GLMs with trend filtering (Eq. (6)) can be optimized using alternating direction method of multipliers (ADMM) (Boyd et al., 2011; Ramdas & Tibshirani, 2016). It can be rewritten as,

$$\begin \min __, ..., \varvec_ } \quad&\sum _1}^ - \ell _(\varvec(g_i) ) + \lambda \Vert D \varvec \Vert _1 \end$$

(9)

$$\begin \iff \quad \min __, ..., \varvec_, } } \quad&\sum _1}^ - \ell _(\varvec(g_i) ) + \lambda \Vert } \Vert _1 \end$$

(10)

$$\begin \text \quad&} - D \varvec = 0 \end$$

(11)

where \(\varvec = (\varvec_^T,..., \varvec_^T )^T\), \(\ell _(\varvec(g_i) )\) is defined in Eq. (5), D represents the difference operator between blocks of \(\varvec\), each block has dimension \(d \times d\).

The augmented Lagrangian is,

$$\begin \begin L_\rho (\varvec, }, })&= \sum _1}^ - \ell _(\varvec(g_i) ) + \lambda \Vert } \Vert _1 + \frac \left\| } - D\varvec + } \right\| ^2 - \frac \Vert } \Vert ^2 \\&= \sum _1}^ - \ell _(\varvec(g_i) ) + \lambda \Vert } \Vert _1 + \frac \left\| } - \sum _1}^ D_ \varvec(g_i) + } \right\| ^2 - \frac \Vert } \Vert ^2 \end \end$$

\(}\) is the scaled dual variable (scaled by \(1/\rho\)). \(\varvec(g_i) \in \mathbb ^d\), \(D \in \mathbb ^\), \(D_ \in \mathbb ^\), \(}, } \in \mathbb ^\). The augmented term is introduced to increase the robustness of the calculation by changing the target into a strict convex problem. Note that \(\rho = 0\) is equivalent to the standard Lagrangian problem. The ADMM update rules are,

Broadcast

$$\begin \varvec(g_i)^ &= \underset(g_i)} \quad -\ell _i(\varvec(g_i) ) + \frac \left\| }^\right. \\ & \quad\left. - \sum _ } D_\varvec(g_j)^\right. \left. - D_ \varvec(g_i) + }^ \right\| ^2, \end$$

(12)

for \(i = 1,...,B\).

Gather

$$\begin }^&= \underset}} \quad \lambda \Vert } \Vert _1 + \frac \left\| } - \sum _1}^ D_\varvec(g_i)^ + }^ \right\| ^2 \end$$

(13)

$$\begin }^&= }^ + }^ - \sum _^ D_\varvec(g_i)^ \end$$

(14)

Equation (12) can be calculated using Newton’s method. Define the target Eq. (12) as \(R( \varvec(g_i) )\), \(\mu _ = \frac \varvec(g_i) \} }\). The gradient and the Hessian matrix are the following,

$$\begin \nabla R &= X_^T (\mu _ - Y_ ) + \rho D_^T \left( \sum _ } D_ \varvec(g_j)^ \right. \\ & \quad\left.+D_\varvec(g_i) - }^ - }^ \right) \end$$

(15)

$$\begin \nabla ^2 R = X_^T \text \left( \mu _ \odot (1-\mu _) \right) X_ + \rho D_^T D_ \end$$

(16)

Equation (13) is equivalent to,

$$\begin }^ = S_\left( \sum _1}^ D_ \varvec(g_i)^1)} - }^ \right) \end$$

where \(S_(\cdot )\) is the coordinate-wise soft-thresholding operator with threshold \(\lambda /\rho\). For j-th entry

$$\begin[S_(\textbf)]_j = \textbf_j - t, & \textbf_j > t \\ 0, & -t \le \textbf_j \le t \\ \textbf_j + t, & \textbf_j < -t \end\right. } \end$$

There are other ways to update the equations above in practice. As suggested by (Boyd et al., 2011, sect. 3.4.5), the algorithm updates each \(\varvec(g_i)\) in turn multiple times before performing the dual variable update.

1.2.2 B.2. Stopping rules

We determine the convergence of the algorithm using primal residuals and dual residuals (Boyd et al., 2011), which stem from the primal feasibility and dual feasibility.

Primal feasibility

$$\begin }^ - D \varvec^ = 0 \end$$

Dual feasibility

$$\begin 0&\in \partial \sum _1}^-\ell _i(\varvec(g_i)^) - \rho D^T (}^/\rho ), \quad } := }^/\rho \\ 0&\in \partial \Vert }^ \Vert _1 + \rho (}^/\rho ), \quad }^ := }^/\rho \end$$

where \(\partial f\) is the subgradient operator of a function f. Note that we use the rescaled ADMM, \(}\) is the original dual variable.

Primal residual

$$\begin }^1)} := }^1)} - D\varvec^1)} \end$$

(17)

Here \(\varvec\) is a stack of \(\varvec(g_i)\).

Dual residual

Since \(}^\) achieves the minimum value of Eq. (13), so

$$\begin 0 &\in\partial \lambda \Vert }^1)} \Vert _1 + \rho \left( }^1)} -D \varvec^1)} + }^ \right) \\ &=\partial \lambda \Vert }^1)} \Vert _1 + \rho }^1)} \end$$

We can see that \(}^\) and \(}^\) always satisfy this part of the dual feasibility. This is also the reason why we set the learning rate as \(\rho\).

As \(\varvec^\) achieves the minimum value of Eq. (12), so

$$\begin 0 \in&\nabla _} \sum _1}^-\ell _(\varvec(g_i)^1)}) - \rho D^T\left( }^ - D \varvec^1)} + }^ \right) \\ =&\nabla _} \sum _1}^-\ell _(\varvec(g_i)^1)}) - \rho D^T\left( }^1)} - D \varvec^1)} + }^ \right) + \rho D^T\left( }^1)} - }^ \right) \\ =&\nabla _} \sum _1}^-\ell _(\varvec(g_i)^1)}) - \rho D^T}^1)} + \rho D^T\left( }^1)} - }^ \right) \\ \Longrightarrow \quad&\rho D^T\left( }^ - }^1)} \right) \in \nabla _} \sum _1}^-\ell _i(\varvec(g_i)^1)}) - \rho D^T}^1)} \end$$

This means that, the following can be viewed as the dual residual.

$$\begin }^ := \rho D^T\left( }^ - }^ \right) \end$$

(18)

1.2.3 B.3. Warm start

ADMM is notorious for slow convergence, especially when \(\lambda\) and \(\rho\) is large. When \(\lambda = \lambda _\), we know the lasso penalty term \(\Vert D \varvec \Vert _1 = 0\) as it shrinks all entries toward zero. So at \(\lambda = \lambda _\) we have,

$$\begin \varvec(g_1) = ... = \varvec(g_B) = \varvec_g^\star \end$$

(19)

All blocks of \(\varvec(g_i)\) are unified. And it achieves the minimum value of the target Eq. (9).

$$\begin \begin&\varvec_g^\star = \underset_g } \sum _1}^ - \ell _(\varvec_g ) \\ \Longrightarrow \quad&\frac_g } \sum _1}^ \ell _(\varvec_g ) = 0 \end \end$$

(20)

At the optimal value, by the stationary condition of \(}^\star\) we also have,

$$\begin }^\star = D \varvec^\star = 0 \end$$

Next we can derive the \(}^\star\) using the stationary condition of \(\varvec^\star\).

$$\begin&\varvec_g^\star = \underset(g_i)} \; -\ell _i(\varvec(g_i) ) + \frac \Vert }^\star - \sum _ } D_\varvec_g^\star - D_ \varvec(g_i) + }^\star \Vert ^2 \\ \Longrightarrow \quad&0 = \frac(g_i)}\left( -\ell _i(\varvec(g_i) ) + \frac \Vert }^\star - \sum _ } D_\varvec_g^\star - D_ \varvec(g_i) + }^\star \Vert ^2 \right) \Bigg |_(g_i) \mathop \varvec_g^\star } \\ \Longrightarrow \quad&0 = - \frac(g_i)} \ell _i(\varvec(g_i) ) \Bigg |_(g_i) \mathop\varvec_g^\star } - \rho D_^T \left( }^\star - D \varvec^\star + }^\star \right) \\ \Longrightarrow \quad&0 = - \frac(g_i)} \ell _i(\varvec(g_i) ) \Bigg |_(g_i) \mathop \varvec_g^\star } - \rho D_^T }^\star \end$$

\(\forall i = 1,...,B\). Now we define,

$$\begin } = \begin -\frac(g_1)} \ell _i(\varvec(g_1) ) \\ ...\\ -\frac(g_B)} \ell _i(\varvec(g_B) ) \end = \begin X_^T (\mu _ - Y_ ) \\ ... \\ X_^T (\mu _ - Y_ ) \end \end$$

where \(X_, \mu _, Y_\) are defined the same as Eq. (15), and the gradient of the PP-GLM log-likelihood function is calculated in the same way. From the stationary condition we know,

$$\begin&\rho D^T }^\star = } \\ \Longrightarrow \quad&}^\star = \frac (D D^T)^D } \end$$

We also need to consider the stationary condition of Eq. (13). As

$$\begin&}^\star = \underset} } \; \lambda \Vert } \Vert _1 + \frac \Vert } - D\varvec^\star + }^\star \Vert ^2 \\ \Longrightarrow \quad&}^\star = S_\left( D\varvec^\star - }^\star \right) = 0\\ \Longrightarrow \quad&}^\star = S_\left( -}^\star \right) = 0, \quad \lambda = \lambda _ \end$$

The last equality must hold as the definition of \(\lambda _\) in Eq. (22) guarantees the zero solution. Then we use the optimal solution \(\^\star , }^\star , }^\star \}\) as the initial values for the ADMM when \(\lambda\) is large. When \(\lambda =\lambda _\), it takes only one iteration to converge of course.

\(\rho\) is an optimization parameter instead of a statistical parameter. Under very general conditions, the ADMM algorithm converges to optimum for any fixed value of \(\rho\) (Boyd et al., 2011). In practice, the rate of convergence and the numerical stability can strongly depend on the choice of \(\rho\) (Ramdas & Tibshirani, 2016). Large \(\rho\) values impose a large penalty on violations of primal feasibility in Eq. (12), so the algorithm favors diminishing the primal residual. Conversely, the definition of \(}^\) in Eq. (18) suggests that small \(\rho\) values reduce the dual residual (Boyd et al., 2011). So we adopt an adaptive strategy to balance the primal and dual residuals as the following,

$$\begin \rho = \tau _} \rho , & \text \Vert }^ \Vert _2> \mu \Vert }^ \Vert _2 \\ \frac} } \rho , & \text \Vert }^ \Vert _2 > \mu \Vert }^ \Vert _2 \\ \rho , &\text \end\right. } \end$$

Since we use a rescaled dual variable, we need to change the \(}\) as well to maintain the same dual variable,

$$\begin }^ = \frac} } }^, & \text \Vert }^ \Vert _2> \mu \Vert }^ \Vert _2 \\ \tau _} }^, & \text \Vert }^ \Vert _2 > \mu \Vert }^ \Vert _2 \\ }^, &\text \end\right. } \end$$

\(\lambda _\) can be derived via KKT conditions (Ramdas & Tibshirani, 2016).

$$\begin 0&\in \frac} \sum _1}^ - \ell _(\varvec(g_i) ) + \partial \lambda \Vert D\varvec \Vert _1 \\ \iff \quad&\frac} \sum _1}^ \ell _(\varvec(g_i) ) = \lambda D^T } \end$$

For some \(}\),

$$\begin }_i \in \, & \text (D\varvec)_i > 0\\ \, & \text (D\varvec)_i < 0 \\ \left[ -1, 1\right] , &\text (D\varvec)_i = 0 \end\right. } \end$$

(21)

So that we get

$$\begin \lambda _ = \Vert (D D^T)^ D } \Vert _\infty \end$$

(22)

where

$$\begin }:= \frac} \sum _1}^ \ell _i(\varvec(g_i) \Big |_=\varvec^\star } ) \end$$

the \(\varvec^\star\) is obtained in Eq. (20).

1.3 C. Stimulus reconstruction

We have presented a method that links channel conductance to specific stimulus filter or post-spike history filter features in time-domain. Next, we will provide a frequency-domain method which is an alternative analysis of how ion channel conductance affects stimulus encoding. We study the frequency properties of the spike-decoded stimulus, which is reconstructed by trained PP-GLMs (Pillow et al., 2008). The decoded stimulus presents the information that has been encoded in the spike train. Stimulus reconstruction provides an intuitive method to investigate how ion channel conductances affect stimulus encoding by comparing the reconstructed stimulus to the actual input stimulus. The method follows the steps in (Tripathy et al., 2013). The stimulus is reconstructed using the maximum a posterior (MAP) estimation of the stimulus given a fitted PP-GLM (Pillow et al., 2008; Tripathy et al., 2013), shown as the following,

$$\begin \max _}}P(} | y_; \, \varvec(g_i), \theta ) = \max _}}P(y_ | s; \, \varvec(g_i)) P(}; \, \theta ) \end$$

(23)

where \(}\) is the vector of full stimulus; \(y_\) is the spike train for the neuron with channel conductance factor \(g_i\); \(\varvec(g_i)\) are the coefficients of the PP-GLM in Eq. (1), \(P(y_ | }; \, \varvec(g_i))\) is the likelihood function given in Eq. (5); and \(P(}; \, \theta )\) is the prior of the stimulus with parameters \(\theta\). As described in section 2.1, the stimulus is the white noise convolved with an alpha function as we introduced in section 2.1. The white noise has a Normal distribution \(N(}, \sigma ^2 I)\). The convolution is a linear transform of the white noise, its corresponding convolution matrix is A. So the prior distribution is \(P(}; \,\theta ) = N(}, \sigma ^2 AA^T)\). We assume the noise variance \(\sigma\) and the alpha function is known.

The optimization of this model is the following. The log posterior of for stimulus reconstruction (Eq. (23)) can be written as,

$$\begin \begin&\log P(y_ | s;\ \varvec(g_i)) P(s; \theta ) \\ &\quad=\sum _1}^ y_ ( K_j s + H_j y_ + \varvec^}) - \sum _1}^ \log (1 + \exp \ + \varvec^}) \}) - \frac s^T (\sigma ^2AA^T)^ s + C \\ &\quad=y_^TKs - }^T \log (1 + \exp \ + \varvec^} \})\\& \qquad - \frac s^T (\sigma ^2AA^T)^ s + C \\ \end \end$$

The above log-posterior is a convex function of s, thus the maximum a posterior (MAP) estimation can be done using Newton’s method. The gradient is,

$$\begin \begin \frac \log P(y_ | s;\ \varvec(g_i)) P(s; \theta ) &=K^T y_ - K^T \sigma \big ( Ks + Hy_ \\ &\quad+ \varvec^} \big ) -(\sigma ^2AA^T)^ s \end \end$$

The Hessian matrix is,

$$\begin \begin \frac \log& P(y_ | s;\ \varvec(g_i)) P(s; \theta ) \\ &=- K^T \text \Big ( \sigma \big (Ks + Hy_+ \varvec^}\big ) \Big ) \\ &\quad K -(\sigma ^2AA^T)^ \end \end$$

K, H are the convolution matrices for the stimulus filter and the post-spike history filter. \(K_j\) and \(H_j\) are the j’th row of the matrices. \(\varvec^}\) is the baseline, it belongs to the parameter \(\varvec(g_i))\). C is the constant which is not a function of s. \(\sigma (\cdot )\) is the sigmoid function. The operators \(\sigma (\cdot )\), \(\exp (\cdot )\) and \(\log (\cdot )\) are element wise, the output has the same dimension as the input. With the gradient and the Hessian matrix, one can use gradient descent method or Newton’s method to get the optimal s. The posterior is a convex function of s, thus it is guaranteed to get the globally optimal solution.

The original stimulus and reconstructed stimulus were compared using the spectrum coherence in different frequency bands. The spectrum analysis was implemented with Welch’s method (MathWorks, 2020). First, the signal was split into overlapping segments. The window length was 256 data points (256 ms width), with 32 data points overlapping. Each window was masked with a Bartlett window. Second, the periodogram was calculated for each window using the discrete Fourier transform, then computed the squared magnitude of the output. All the periodograms were then averaged. Third, we estimated the magnitude-squared coherence, which is a function of frequency with values between 0 and 1, indicating how well the input signals x matched to y at each frequency. The estimator for the coherence is the following (Kramer, 2013),

$$\begin \hat_(f)=_(f)|^}_(f)\hat_(f)}} \end$$

(24)

where \(\hat_(f)\) is the estimated cross-spectral density between x and y, \(\hat_(f)\) and \(\hat_(f)\) are the estimated auto-spectral density. The estimated spectral densities were estimated by averaging the periodograms of all windows.

Fig. 8

Stimulus reconstructions and spectral coherence. A An example stimulus reconstruction for conductance scaling of 1.5, 1.0, and 0.5 (colored lines) compared to the actual stimuli (gray line) for the MC KAchannel. B Magnitude squared coherence between the stimulus reconstruction and the mean stimulus for conductance scaling of 1.5, 1.0, and 0.5. C The difference in coherence between the conductance scaling and control scaling of 1.0. D-F The mean coherence across indicated frequency bands as a function of conductance scaling factor. The gray dotted line represents control scaling factor. All panels are for the MC KAchannel

Another way to examine the PP-GLM, which is an encoding model, is through decoding. Decoding is the process of estimating a reconstruction of the original stimulus given a spike train and a trained PP-GLM (Eq. (23); Fig. 8A). We then compare the reconstructed stimulus to the original stimulus by measuring the coherence as a function of the signal frequency (Fig. 8B, C). We consider only stimulus reconstructions from PP-GLMs trained with optimal trend filtering penalty hyperparameter, as the reconstructed stimuli for PP-GLMs trained without trend filtering were nearly identical (data not shown). This is expected, as the goodness-of-fit is nearly identical between \(\lambda = 0\) and \(\lambda = \lambda ^*\) (Fig. 5C). The coherence analysis allows estimation of specific frequency components that are, or are not, encoded when scaling different ion channel conductances (Fig. 8B). Here we evaluate how ion channel conductance scaling affects the coherence between the reconstructed stimulus and the original stimulus, by measuring the difference between scaled ion channel conductances and the control ion channel conductance (Fig. 8C). For example, when scaling the MC KAchannel, increasing KAchannel conductance generally reduces coherence across the frequency spectrum, whereas decreasing MC KAchannel conductance shows increased coherence at specific frequencies 35-50 Hz and 70 Hz (Fig. 8C). Generally, the coherence measures are fairly noisy, which we can smooth by averaging over well characterized frequency bands (Fig. 8D-F). The MC KAchannel conductance scaling affects the encoding of mid range, beta frequencies (Fig. 8D-F), with only moderate effects on low range, theta frequencies and high range, gamma frequencies. This suggests a prominent role for the MC KAchannel in the encoding of mid range, beta frequencies. Overall, the additional approach of examining stimulus reconstructions further reveals how different ion channel conductance scaling affects the encoding of specific stimulus features.

1.4 D. Inverse mapping

The main text focuses on the mapping from the biophysical model parameter space to the statistical model parameter space. This section discusses the inverse mapping: how to predict biophysical properties using observed spike trains. Let \(\varvec_}\) be the parameters of the biophysical model and \(\varvec_}\) be the parameters of PP-GLM, which can be estimated with some level of uncertainty. A natural method of implementing the inverse mapping is building the probability \(p(\varvec_} | \varvec_})\), where \(\varvec_}\) comes from the GLM fitted to spike trains with unknown biophysical properties.

Before inferring the biophysical parameter, \(p(\varvec_} | \varvec_})\) needs to be learned, which can be approximated as,

$$\begin \begin&p(\varvec_} | \varvec_}) \propto p(\varvec_}, \varvec_}) = p(\varvec_} | \varvec_}) p(\varvec_}) \\ &\approx\sum _ p(\varvec_} | \textrm_l^}, \varvec_}) p(\textrm_l^} | \varvec_}) p(\varvec_}) \end \end$$

(25)

where \(\textrm_l^}\) are samples indexed by l coming from the biophysical simulator with biophysical parameter \(\varvec_}\). This requires a biophysical simulator for model training. Equivalent real experiments need to collect a training dataset with both spike trains and the underlying biophysical properties. In the above equation, marginalizing out the spike trains is approximated by summing over generated spike trains. The prior \(p(\varvec_})\) is assumed to be a uniform distribution within a reasonable range. Given spike trains and biophysical parameters, corresponding GLM parameters can be calculated using the method in the main paper. In summary, we can draw samples \((\varvec_}, \varvec_})\) from their joint distribution with a certain approximation. But this can not be directly used to infer \(\varvec_}\) given arbitrary \(\varvec_}\), that needs one more step to approximate the probability through some parametric form. Here is one example,

$$\begin p(\varvec_} | \varvec_}) \propto p(\varvec_}, \varvec_}) \propto \exp \left\_} - \theta ^T \varvec_} \Vert _F^2 \right\} \end$$

(26)

where \(\Vert \cdot \Vert _F\) is the Frobenius norm. Approximating the above distribution is equivalent to training \(\theta\) with the samples drawn from \(p(\varvec_}, \varvec_})\).

Finally, we run a simple experiment to estimate MC KAconductance with the linear approximation. The estimation error of the channel achieves 0.01 (the unit is scaling factor).

1.5 E. Goodness-of-fit Fig. 9

Goodness-of-fit test. The dashed line is the 95% uniform CI. The figure shows an example of MC KAchannel. Different subplots corresponds to different scaling factors indicated by the title. The curve staying within the CI implies a good fit

The goodness-of-fit of the model is evaluated using the KS test based on the time-rescaling theorem (Kass et al., (2014; Haslinger et al., 2010; Brown et al., 2002). The results of MC cell KAchannel is shown in Fig. 9. The model shows good fit for all different scaling factors.

Fig. 10

A comparison between biophysical data of MC cell KAchannel in panel A and PP-GLM simulated spikes in panel B

Besides the KS test, we also compare the biophysical spike trains and the spike trains generated by the fitted PP-GLM as shown in Fig. 10. Both the spike train raster plots and the PSTHs show a good match.