In this study, we used the experimental data published in (Lei et al., 2019, https://github.com/CardiacModelling/hERGRapidCharacterisation). Therefore, the hERG model and basic settings we use are the same as those of (Lei et al., 2019). For ease of training and evaluation of our deep learning model, we excluded 11 of the 211 cells, which seemed to have a large difference between the experimental results and the results produced by the hERG model. The currents were recorded for the “staircase protocol” (Figure 2). As shown in the top image in Figure 2, each step is 500 ms, long enough to see the characteristics of I K r . Thus, it is possible to observe the dynamics at different voltage values. Lei et al. showed that their protocol provided enough information to infer true parameters through a synthetic data study (Lei et al., 2019).

This hERG model with Hodgkin and Huxley style structure used in this experiment is Beattie’s model (Beattie et al., 2018) slightly improved by ten Tusscher et al. (ten Tusscher et al., 2004), and I _Kr is the same as Eq. (1). I K r = g K r a r V − E K (1) where g K r is the maximum conductance, and a and r denote a Hodgkin and Huxley activation gate and an inactivation gate, respectively. V is the transmembrane voltage. E K is called Nernst potential or the reversal potential and obtained by Eq. (2). E K = R T z F ln K + o K + i (2) where R   i d e a l   g a s   c o n s t a n t : 8.314472   J K ∙ m o l T   a b s o l u t e   t e m p e r a t u r e : 298.15   K F   F a r a d a y ′ s   c o n s t a n t : 96485.3415   C m o l z   v a l e n c y   o f   t h e   i o n s : + 1   f o r   K + K + o   e x t r a c e l l u l a r   c o n c e n t r a t i o n : 4   m M K + i     i n t r a c e l l u l a r   c o n c e n t r a t i o n : 110   m M

The model has nine parameters θ = g K r , p 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 , p 8 , and Figure 3 shows its structure, where d a d t = a ∞ − a τ a , d r d t = r ∞ − r τ r a ∞ = k 1 k 1 + k 2 , r ∞ = k 4 k 3 + k 4 τ a = 1 k 1 + k 2 , τ r = 1 k 3 + k 4 k 1 = p 1 ⁡ exp p 2 V k 2 = p 3 ⁡ exp − p 4 V k 3 = p 5 ⁡ exp p 6 V k 4 = p 7 ⁡ exp − p 8 V where k 1 is an activation rate, k 2 is a deactivation rate, k 3 is an inactivation rate, and k 4 is a recovery rate.

A lot of data is required to train. The best condition for good performance is the presence of a large amount of experimental data. However, there are only 200, and it is very insufficient for learning and testing with them. Even if a large amount of data exists, fitting work for labeling requires a lot of time and resources, which is contrary to the purpose of our study. So, we generated a large amount of data using simulations of the hERG channel to compensate for the lack of data. The dataset consists of I _Kr as input data and nine parameters θ of the hERG model as the target data. Five hundred thousand data generated by simulation were used for training. For 200 experimental data, after sorting by name, odd-numbered data were configured as validation dataset and even-numbered data were configured as test dataset.

First, parameters for the hERG channel were generated under the following two conditions. The first condition is that each parameter follows a uniform distribution within a specific range, as in Eq. (3). g K r ∼ U g min , g max p 1 , p 2 , p 3 , p 4   ∼ U a , b (3) p 5 , p 6 , p 7 , p 8   ∼ U c , d where U ∙ represents a uniform distribution. In this study, for the conductance g K r , g min = 100 [pA/V] and g min = 500000 [pA/V]. For p 1 , p 3 , p 5 , and p 7 , a and b are 0.0001 [s⁻¹] and 10⁶ [s⁻¹], respectively. For p 2 , p 4 , p 6 , and p 8 , c and d are 0.0001 [V⁻¹] and 400 [V⁻¹], respectively. The second condition is that each parameter must satisfy the inequalities of Eq. (4): 0.0167 < p 1 exp p 2 * V max < 10 6 0.0167 < p 3 exp − p 4 * V min < 10 6 (4) 0.0167 < p 5 exp p 6 * V max < 10 6 0.0167 < p 7 exp − p 8 * V min < 10 6 where V min = − 0.12 and V max = 0.06 . The lower and upper bounds of the above conditions were determined by the constraints of physical and physiological phenomena (Beattie et al., 2018). Figure 4A shows the distribution of each parameter for 100,000 data generated by the above conditions.

Second, we used simulations to generate the current I K r corresponding to the parameters sampled above. Myokit (Clerx et al., 2016) with CVODE solver (Hindmarsh et al., 2005) was used for the simulation. The tolerance settings for CVODE were abs_tol = 10^–8 and rel_tol = 10^–10, as in the condition in (Lei et al., 2019). The length of the experimental data was 15.4 s with a sampling rate of 5 kHz. In this study, we reduced the data number to 1/50 with sampling one point every 100 because there seemed to be no problem reflecting the trend. Because we only used the data generated by simulations for training, the results of this study depend on the similarity between the simulation data and experimental data. Therefore, noises, α , were added to the simulation data. The noises were extracted from a normal distribution α ∼ N 0 , σ 2 . σ is 10.84, which was measured at the steady-state current in the experimental data. I K r exp ≈ I K r s i m + α (5) where I K r e x p and I K r s i m represent the experimental and simulated current, respectively.

Our method involves several simple preprocessing processes on the data for learning. In Figure 4A, The distributions of p 1 , p 3 , p 5 and p 7 are very clustered in a specific range. To make them as uniform as possible, the log scale was applied. The min–max normalization was then used to transform all parameter ranges between 0 and 1, as shown in Figure 4B.

Recurrent neural network (Sherstinsky, 2020) and LSTM (Hochreiter and Schmidhuber, 1997; Sherstinsky, 2020) series models have been widely used to analyze time-series data, such as the current data that we want to analyze, and recently, transformer-based models (Vaswani et al., 2017; Lin et al., 2022) are leading this field. Since CNN is designed for the purpose of extracting information between adjacent values of data, it obtains spatial information well in the local domain (Krizhevsky et al., 2017). Transformer calculates the relationship between all elements of input data through attention, so it understands overall features better than CNN, but is weaker than CNN in extracting local information, and requires a very large size dataset for this purpose (Dosovitskiy et al., 2020). We thought our goal was closer to finding changes in the characteristics and patterns at specific times than predicting the current change over time. Therefore, we adopted a CNN-based model rather than a transformer-based model, and among them, an EfficientNet (Tan and Le, 2019) type model. We called our model hERGNet. Our hERGNet is very simple, as shown in Figure 5. It consists of a CNN-encoder network and a pretrained EfficientNetV2-M (Tan and Le, 2021). The CNN-encoder network extracts the features of the spectrogram and increases the number of channels to three to obtain an image-like shape, thereby making it possible to use the pretrained EfficientNetV2-M. EfficientNetV2-M is responsible for finding parameters by extracting features from encoded data.

We converted the current into a spectrogram with a frequency perspective, as shown in Figure 6, so that hERGNet can better learn the characteristics of the current data. By adding frequency features to current data consisting of only time and intensity, it has the advantage of increasing information about data and transforming it into a two-dimensional form like an image, making it easier for a 2D CNN-based model to learn. In this study, the parameters “n_fft,” “hop_length,” and “win_length” of Short-Time Fourier Transform (Owens and Murphy, 1988) for spectrogram transformation were set to 256, 12, and 48, respectively.

Mean Squared Error (MSE) was used as the loss function, so the cost function J(θ) is the same as Eq. (6) J θ = 1 N ∑ i = 1 N f θ x i − y i (6) where x i is the spectrogram of the current I K r , f θ is the neural network, and y i is the true parameters of the ion channel.

In general, it is known that the higher the resolution of the input, the higher the performance. However, the higher the resolution, the more memory is required, which is time-consuming to train due to the small batch size. Also, if you train all the generated 500,000 data from the beginning, the time will increase even more. Therefore, we first trained the hERGNet with a small resolution and a small number of data and then proceeded with transfer learning by increasing the resolution and the number of data. First, learning was performed on 300,000 spectrograms with a 97 × 97 resolution to 200 epochs. Then, we increased the resolution to 129 × 129 and performed transfer learning to 140 epochs using all 500,000 data.

All fitting operations were the same as those in (Lei et al., 2019), and their open-source library PINTS (Clerx et al., 2019) was used. Furthermore, the CMA-ES algorithm (Hansen, 2006; Khan, 2018) was used as a global optimization algorithm to fit the model to the experimental data, and Markov Chain Monte Carlo (Jasra et al., 2007) with adaptive Metropolis (Haario et al., 2001) was used to explore the posterior probability distribution. Parameter inference was performed three times based on the initial value for each cell. The first was when the initial value was given as a parameter predicted by our hERGNet, the second was a prior value, and the last was given as a value randomly extracted from the previously described parameter distribution.