A Pairwise Connected Tensor Network Representation of Path Integrals

It has been recently shown how the tensorial nature of real-time path integrals involving the Feynman-Vernon influence functional can be utilized using matrix product states, taking advantage of the finite length of the non-Markovian memory. Tensor networks promise to provide a new, unified language to express the structure of path integral. Here, a generalized tensor network is derived and implemented specifically incorporating the pairwise interaction structure of the influence functional, allowing for a compact representation and efficient evaluation. This pairwise connected tensor network path integral (PCTNPI) is illustrated through applications to typical spin-boson problems and explorations of the differences caused by the exact form of the spectral density. The storage requirements and performance are compared with iterative quasi-adiabatic propagator path integral and iterative blip-summed path integral. Finally, the viability of using PCTNPI for simulating multistate problems is demonstrated taking advantage of the compressed representation.

I. INTRODUCTION
Tensor networks (TN) are designed to be compact "factorized" representations of high-ranked tensors. Probably the most common use of TN in physics is related to representations of the quantum many-body wavefunction which, in general, is also a high-ranked tensor. This use has been widely demonstrated in a multitude of methods such as the density matrix renormalization group (DMRG) [1,2] which uses a Matrix Product State (MPS) [3,4] representation, and multiconfiguration time-dependent Hartree (MCTDH) [5] and its multi-layer version (ML-MCTDH) [6][7][8] which use tree tensor networks. For multidimensional systems, an "extension" of MPS to multiple dimensions called projected entanglement pair states (PEPS) [9] is used. For systems at critical points, an MPS representation does not work because of long-range correlations necessitating the use of the so-called multi-scale entanglement renormalization ansatz (MERA) [10,11]. Tensor networks, since its introduction, have proliferated in various diverse fields requiring the use of compact representations of multidimensional data like machine learning and deep neural networks.
While quantum dynamics at zero temperature can often be simulated using wave-function based methods like time-dependent DMRG [2,12,13] or MCTDH, at finite temperatures, owing to the involvement of a manifold of vibrational and low frequency ro-translational states in the dynamics, they suffer from an exponentially growing computational requirements. Feynman's path integral provides a very convenient alternative for simulating the time-dependent reduced density matrix (RDM) for the system. The vibrational states of the "solvent" introduced as harmonic phonon modes under linear response [14] are integrated out leading to the Feynman-Vernon influence functional [15]. Identical influence functional also arises in dealing with light-matter interaction through the integration of the photonic field.
The primary challenge in using influence functionals and path integrals is the presence of the non-local history-dependent memory that leads to an exponential growth of system paths. While many recent developments have helped improve the efficiency of simulations [16][17][18][19][20], each of them utilize very different and deep insights into the structure of path integrals. It has recently been shown that the MPS representation can be very effectively utilized to reformulate real-time path integrals involving the influence functional leveraging the finite nature of the non-local memory [21][22][23][24]. While the MPS structure is the simplest tensor network that can be used, the 1D topology is probably not optimal when the non-Markovian memory spans a large number of time-steps and suffers from growing bond dimensions. In this paper, an alternate generalized tensor network that directly captures the pairwise interaction structure of the Feynman-Vernon influence functional, is introduced. This pairwise connected tensor network path integral (PCTNPI) has an extremely compact representation, that can be efficiently evaluated, allowing us to go to much longer non-Markovian memories without resorting to various techniques of path filtration. Tensor networks show great promise in being a unifying language for formulating and thinking about path integral methods.
The construction and evaluation of the tensor network is discussed in Sec. II. In Sec. III, we illustrate some typical applications of the algorithm. The memory usage is also reported for various parameters. The implementation of this method utilized the open-source ITensor [25] library for tensor contractions allowing for extremely efficient tensor contractions using highly efficient BLAS and LAPACK libraries. We end the paper in Sec. IV with some concluding remarks and outlook on future explorations.

II. METHODOLOGY
Consider a quantum system coupled to a dissipative environment described by a Caldeira-Leggett model [26][27][28] whereĤ 0 is the Hamiltonian of the D-dimensional system of interest shifted along the adiabatic path [29]. If the quantum system can be described by a two-level Hamiltonian, thenĤ 0 = σ z − Ωσ x , whereσ z andσ x are the Pauli matrices.Ĥ env represents the Hamiltonian of the reservoir or environment modes which are coupled to some system operatorŝ. The strength of the j th oscillator is c j . While we are using a time independent Hamiltonian for simplicity, time-dependence from an external field in the system Hamiltonian can be captured through the corresponding system propagator in a straightforward manner. For a problem where the environment is in thermal equilibrium at an inverse temperature β = 1 k B T , and its final states are traced out, the interactions between the system and the environment is characterized by the spectral density [14,26] In fact, the spectral function, S(ω) corresponding to the collective bath operator X = − j c j x j is related to the spectral density as follows [30]: For environments defined by atomic force fields or ab initio calculations, it is often possible to evaluate the spectral density from classical trajectory simulations. The dynamics of the RDM of the system after N time steps, if the initial state is a direct product of the system RDM and the bath thermal density is given as: where Here,Û is the short-time system propagator for ∆t and {s + j } and {s − j } are the forward-backward system paths. The Feynman-Vernon influence functional [15], F [s ± j ], is dependent on the system path s ± j and the bath response function that is discretized as the η kkcoefficients [31,32]. The influence functional depends upon the history of the system path, leading to the wellknown non-Markovian nature of system-environment decomposed quantum dynamics. Notice that it can be factorized based on the "range" of interaction in the following manner: The influence functional creates pairwise interactions between points that are temporally separated. As it has been shown, if MPS and MPO are used to model the influence functional, the fact that these interactions can spread across long temporal spans leads to an increase in the effective bond dimension. Here, the goal is to create a structure that naturally and efficiently accounts for the pairwise interactions that span long temporal separations while not being associated with any one particular representation.
To motivate the tensor network representation, first consider the Markovian part of Eq. (5), involving just the propagators and the terms of the influence functional coupling consecutive time points. These terms can be simply rearranged as: Here, we are implicitly summing over repeated indices that do not appear on both sides of the equation. The ± labels on the site indices of the tensors are omitted for convenience of notation. The superscript, 1, on P is there to denote the maximum distance of interaction that we have incorporated. Equation (9) is already a tensor network; more specifically it is series of matrix multiplication as shown in Fig. 1. Let us now bring the "next-nearest neighbor" interactions I (2) sj−2,sj . Clearly, it is not possible to directly contract the I (2) sj−2,sj tensor to the P (1) tensor because the internal s j 's have already been traced over. To make it possible to incorporate the I (2) tensors, we augment the K tensors as follows: It is convenient to think of lower indices as the "input" indices and the upper indices as the "output" indices, though there is no other mathematical significance to the positioning of the indices. With this input-output convention in mind, it is easy to see that the internal augmented K tensors duplicate and flip the order of the input indices, s j . This ensures that indices that differ by two time steps are now placed adjacent in the output layer. Now, the Markovian terms and the I (2) terms can be combined and we get: which is depicted in Fig. 2. Notice that the index s j still connectsK sj−1,sj andK sj ,sj+1 , as in P (1) , but now there is another connection that goes through the I (2) tensor in a "triangular" form. This feature of an I (α) with a higher α acting as a bridge betweenK or I tensors with smaller values of α would become a recurring motif in this tensor network. The pattern for inclusion of the rest of the non-local interactions is quite similar. Note that in Fig. 2, if we did the same "trick" of duplicating and flipping the order of the inputs, in the next layer indices that differ by three time points, like s 0 and s 3 , s 1 and s 4 , are going to be adjacent. Hence, this can now be multiplied sj−3,sj . Continuing like this, we can complete the network. The diagram is shown in Fig. 3 [33], these augmented tensors are going to be written as K, I (2) , . . . . The tensor network shown in Fig. 3, which we will conventionally denote by P (∞) , represents the final Green's function for the propagation of the system RDM having incorporated the non-local influence from the environment. So, ρ(t) = P (∞) ρ(0).
If the system is defined to have D states, then in Fig. 3, all the indices have D 2 dimensionality corresponding to each of the possible combination of forward-backward states. However, this is not optimal. Notice that the influence functional tensors, I (α) s k−α ,s k for a time difference of α, depends only on the "difference" coordinate, ∆s k = s + k − s − k of the latter time point. So, currently, we are carrying over more information than we need to.
To take care of this redundancy, we need to redefine the K tensors to not just duplicate the input indices, but to project the "latter" index onto its difference coordinates as follows: Notice that the upper-right output indices in the diagrams remain exactly the same. Only the upper-left output index of K changes. Therefore, the tensor K r1 s0,s1 remains unchanged. The dimensionality of the "l" indices is the number of unique values of ∆s = s + − s − that the system can have. For a general D-level system, this value is B = D 2 − D + 1 instead of D 2 , however the actual symmetries present in the system might reduce this even further. Finally, the influence functional tensors have to be changed to be consistent, viz.Ĩ (α) these changes, the basic topology of the network remains the same. The new network with the different dimensions is shown in Fig. 4.
Having discussed the tensor network, now let us turn to the job of contracting it. Typically, many tensor networks are constructed using singular value decomposition (SVD) and evaluated via the truncation of the singular values [21,24]. The PCTNPI network is constructed without resorting to any SVD calculations and consequently "exact." The goal now is to find an optimal contraction scheme that preserves this "exactness." The storage cost, S, is also evaluated at the end of every step. The canonical contraction order that we discuss below has been marked out in cyan arrows in Fig. 4. For a simulation with N time steps: s0,s N and contract it withĨ  9. Finally, the last tensor,K s N −1 ,s N is contracted.

Multiply allĨ
In the above contraction scheme, we multiply the initial condition, ρ s0 , and get the final RDM. While this leads to a more efficient algorithm in terms of the storage and computational cost, it is possible to reformulate the scheme in terms of the Green's function by not involving the initial condition in the contractions and evaluating P (∞) . An in-depth analysis of the memory and computational cost is given in Appendix A. Of course, the storage requirement grows to a maximum of (D 2 ) 2 B N −2 before decreasing continuously. This naïve contraction scheme does not solve the problem of storage. Still, as would be illustrated in Sec. III, PCTNPI outperforms both traditional iterative quasi-adiabatic propagator path integrals (QuAPI) [31,32], and iterative blip summed path integral (BSPI) [17,18], when used without path filtration, in the memory lengths that can be accessed without any sort of filtration. In a future work, filtration schemes on top of PCTNPI would be introduced that can not only deal with this problem, but would also avoid the construction and storage of the full tensor. The focus of this paper is however on the tensor network and its performance in the most naïve implementation.
It is well-known that the non-local memory of the influence functional dies away with the distance between the points, allowing for a truncation of memory. This idea is commonly used both in Nakajima-Zwanzig generalized quantum master equations [34][35][36] and iterative QuAPI [31,32]. In the framework of PCTNPI, the length of the non-Markovian memory is equal to the depth of the resultant network. The topmost tensor encodes the interaction between the most distant points, while the bottom most tensor captures the Markovian interactions coming through the propagator and the I (1) terms.
At two time-steps of memory, that is L = 2, we basically get Fig. 2. In Fig. 5, we show the structure of the network for a 5-step propagation with L = 3. Because s 0 does not interact with s 4 or s 5 , it is not necessary to store and evaluate the full diagram at once, but it can be built iteratively. The first edge, corresponding to interactions with s 0 is contracted, and multiplied by the second edge, using the canonical contraction scheme discussed previously. As soon as this is done, the storage of the first edge can be freed, and the third edge can be contracted. This iteration scheme turns out to be identical to the iteration scheme in iterative QuAPI. The first steps of the iteration algorithm is pictorally outlined in Fig. 6. Makri [17] has shown that it is possible to think of the memory as arising from two different causes. The influence functional F can be rewritten in terms of the real and imaginary parts of the η-coefficients as: where ∆s k = s + k − s − k ands k = 1 2 (s + k + s − k ). The part of the influence functional that arises from Re η kk is called the classical decoherence factor. It corresponds to stimulated phonon absorption and emission [37]. This can also be obtained through classical trajectory-simulations and reference propagators [38] in a Markovian manner. All effects of temperature is captured in the classical de-coherence term. The term with the Im η kk is the backreaction that leads to quantum decoherence. This part of the memory is truly non-local and temperature independent.
As a cheap approximation to the dynamics, it is possible to do a simulation with classical decoherence, that would become increasing accurate as the temperature of the simulation rises. In this, the full η kk coefficients are used only when k = k or k = k + 1, and otherwise the imaginary part of η kk is ignored. (Actually, the true expressions for classical decoherence would include the full η kk coefficients only when k = k and the real part otherwise. In PCTNPI, we can include the case of k = k + 1 as well at the same storage and computational cost.) Effectively, we are modifying the I (α) s k ,s k operators to be exp (−1/ Re(η kk )∆s k ∆s k ) when α = k − k ≥ 2. Just like before when the s k lines carried unnecessary information, now the s k lines carry more information than they need to. We only need to know about ∆s k . Thus we can make the required changes to the dimensionality of the indices by putting in the corresponding projector operators in the K tensors, thereby reducing the cost of computation even further. The network for the classical decoherence simulations would have exactly the same structure as Fig. 3 with all edges except the base ones being B dimensional. This approximation is especially accurate at short times.

III. RESULTS
As illustrative examples, we apply PCTNPI to a twolevel systems (TLS) coupled bilinearly to a dissipative environment:Ĥ The dissipative environment is chosen to be defined by Ohmic model spectral densities, which are especially useful in modeling the low frequency ro-translational modes. We use the very common Ohmic form with an exponential cutoff, where ξ is the dimensionless Kondo parameter and ω c is the characteristic cutoff frequency and the Ohmic form with a Drude cutoff, where κ is a measure of the coupling strength. Generally these model spectral densities are often thought to be fully characterized by a reorganization energy (22) and the cutoff frequency, ω c . The reorganization energies for the exponential and the drude cutoff spectral densities are as listed below: As we demonstrate through the examples, though the reorganization energy and the cutoff frequency are same, the exact dynamics of the reduced density matrix is highly dependent on the form of the "decay function." Consider a symmetric TLS ( = 0) and Ω = 1 interacting strongly (ξ = 2) with a sluggish bath (ω c = Ω) initially localized on the populated system state 1. The bath has a reorganization energy of λ = 4 and is held at an inverse temperature of Ωβ = 1. The dynamics was converged at ∆t = 0.125, and a memory length L = 16. The convergence is shown in Fig. 7 (a) for an Ohmic bath with an exponential decay. Full quantum-classical simulations for this parameter is available [39]. If the Drude form of decay is used, the dynamics changes quite signif-icantly. The comparison between the dynamics arising from the two spectral densities is shown in Fig. 7 (b). Next, consider a case where not only is the dynamics different between the two different decay functions, but the converged non-Markovian memory length is different as well. The dynamics of the same TLS as above ( = 0, Ω = 1) is now simulated in a bath with the reorganization energy λ = 8 and a characteristics cutoff frequency ω c = 5. The bath is equilibrated at an inverse temperature of Ωβ = 5. The time-step is converged at Ω∆t = 0.125. The convergence of the dynamics of the reduced density matrix on changing the memory length, L, is shown in Fig. 8. While the memory length for the exponential decay function spectral density is quite close to convergence at L = 14, for the Drude spectral function, it converges at L = 10.
In Fig. 9, we consider a TLS coupled to a strongly coupled Ohmic bath with an exponential cutoff (ξ = 1.2, ω c = 2.5Ω) equilibrated at a high temperature Ωβ = 0.2. The converged time step is ∆t = 0.125. The classical memory calculations converge at a comparatively lower memory length, L and agree quite well with the full simulations at short times. Though at intermediate and long times, the classical decoherence dynamics differs from the true dynamics, this can often be enough for estimating timescales of processes, especially using rate theory [40][41][42][43].
Next, the storage requirements of PCTNPI is compared with that of iterative QuAPI and iterative BSPI in Fig. 10 (a). To keep the comparisons fair, the iQuAPI and iBSPI methods were run without any path filtering. It is quite clear from the plot that the scaling of PCT-NPI is essentially "like" that of iBSPI, i.e. for a TLS, 3 L scaling for iBSPI and PCTNPI vs 4 L scaling of iQuAPI. However, the prefactor is much smaller, allowing us to access much longer memories with limited resources. In fact, this difference in the prefactor would grow with the dimensionality of the quantum system. In iBSPI, there would be 3 L paths for a memory length of L, but the storage is more than just a number corresponding to each path. It stores a small dimensional matrix for each path. This is the cause of the larger prefactor.
A comparison of the run times of PCTNPI with respect to iBSPI without any filtration for a simulation of 100 time steps is presented in Fig. 10 (b). A laptop with Intel® Core TM i5-4200U CPU with a clock speed of 1.60GHz was used for these benchmark calculations. These measurements are not going to be consistent with similar benchmarks run on other machines, but the basic trends would continue to hold. The PCT-NPI algorithm is built on top of ITensor and automatically uses parallel BLAS and LAPACK wherever possible. There is no standard iBSPI code. The iBSPI program used for these benchmarks was manually parallelized with OpenMP loop parallelization.
As a final example, consider a molecular wire described by the tight-binding Hamiltonian involving N sites:   The site energy of the j th site is j and the nearest neighbor couplings are V . The sites are separated by unit distance such that |σ j are eigenstates of the position operator,ŝ |σ j = (j − 1) |σ j . The site energy of all but the first site is chosen to be zero j = 0 for j = 1 and 1 = 1. The intersite coupling is chosen to be V = 0.025 [44]. The computational cost grows exponentially with the number of sites. To test the efficiency of the basic contraction scheme outlined here, we use a system with N = 4 sites. The bath is characterized by an Ohmic spectral density with an exponential cutoff, Eq. (20) with ω c = 4 and ξ = 0.12 [24] equilibrated at an inverse temperature of β = 0.1. As discussed in Sec. II, the scaling of the algorithm would go as B < D 2 . The symmetry of the Hamiltonian in this case ensures that the number of unique values of ∆s, B = 7 for this 4 state system, which is even less than the D 2 − D + 1 for a completely general Hamiltonian. The population dynamics of all the states is shown in Fig. 11. An initial state with only the first site populated was used. Because of the high temper- ature of the bath, the classical decoherence simulation produces practically identical dynamics but converges at a smaller memory length L.

IV. CONCLUSION
A novel tensor network is introduced to perform path integral calculations involving the Feynman-Vernon influence functional. This pairwise connected tensor network path integral (PCTNPI) captures the pairwise interaction structure of influence functional. PCTNPI can be contracted efficiently, and minimizes the storage requirements as far as possible without resorting to various path filtration algorithms. Iterative decomposition of the memory is also possible in an elegant manner. Comparisons between PCTNPI and iQuAPI and iBSPI show the scaling of memory requirements of PCTNPI to be similar to iBSPI, but much smaller.
PCTNPI provides an alternative to the MPS representation [21,24], serving as a small step in further elucidating the deep relation between tensor networks and path integrals. While no path filtration scheme has been developed, PCTNPI is already quite useable. It can easily incorporate classical trajectories through harmonic backreaction quantum-classical path integrals [37] thereby making it possible to include anharmonic effects of the environment in an approximate manner without any additional cost. Additionally, harmonic backreaction also leads to an increase in the converged time-step and a decrease in the effective memory length such that some ultrafast reactions can be simulated directly. Taking advantage of the extended memories that are accessible with PCTNPI, the combined method would be able to simulate systems with strongly coupled sluggish realistic solvents with high reorganization energy. This promises to be a fruitful avenue of research in terms of applications (a) Contraction of internal tensor along the left edge (j ≥ 2). to electron and proton transfer reactions.
Algorithms based on MPS representations of the augmented reduced density tensor [21] or of the pathdependent Green's function [24] can be thought of as particular optimized re-factorizations of the PCTNPI network. We have demonstrated the viability of evaluating PCTNPI in a brute force manner compared to other methods. This suggests that using the PCTNPI network directly to generate other optimized representations might also lead to novel methods.
While ideas of path filtration were not a consideration of the present paper, schemes based on the singular value decomposition (SVD) can be incorporated with PCTNPI, leading to a method that significantly reduces the storage, since the full tensor would not need to be computed and stored. This development would be discussed in a future publication. cost of contraction, the left "edge" of the triangular network is first considered. Consider contractingĨ (j) s0,sj , for j ≥ 2, with two D 2 indices and one B index, as schematically indicated in Fig. 12 (a). The part that has already been contracted has one D 2 index and (N − j) B indices. Therefore, the cost of contraction is (D 2 ) 2 B N −j+1 . The space requirement at this stage is D 2 B N −j+1 . To finish the contraction of the left-most edge of the triangle, we need to multiply byK s0,s1 ρ s0 leading to the tensor network shown in Fig. 12 (b). The resultant tensor does not have a index corresponding to s 0 because that has been traced over. The computational cost of this step is C = (D 2 ) 2 B N −1 and the storage becomes S = D 2 B N −1 . Now, the second parallel edge is to be contracted. This step however is started from the bottom, i.e. fromK s1,s2 .
The first contraction, shown in Fig. 13 (a), is the most costly step in the entire algorithm. The computational cost of this step is C = (D 2 ) 3 B N −1 and the storage requirement increases to S = (D 2 ) 2 B N −2 . Continuing with the other intermediate tensors of the first parallel edge, notice that the cost of contraction remains constant at C = (D 2 ) 3 B N −1 and the space required remains constant at S = (D 2 ) 2 B N −2 . Finally, the last, top-most tensor of this edge is to be contracted. This is illustrated in Fig. 13 (b). The computational cost is C = (D 2 ) 2 B N −1 . The storage cost now drops to S = D 2 B N −2 . Now, consider contracting a general diagonal edge, say the j th one. The resultant tensor from the previous contraction has one D 2 index and N − j + 1 B indices. Contracting theK tensor leads to a tensor with two D 2 indices and (N − j) B indices. The cost of this contraction is C = (D 2 ) 3 B N −j+1 and the storage is S = (D 2 ) 2 B N −j . For all the intermediate tensors at this stage, once again both the computational costs and the storage costs remain the same. On contracting the last tensor of this diagonal, the storage drops to S = D 2 B N −j .
Below we list the total computational cost for contracting each of the "parallel" edges. The edge number is given as the subscript.
The prefactor of the computational and storage costs is lower for classical decoherence simulations: It goes from a power of D 2 to the corresponding power of B. It is clear that the complexity of the entire contraction goes as O B N −1 and the peak storage requirement is O B N −2 .


Time-frequency time-space long short-term memory networks for image classification of histopathological tissue

Image analysis in histopathology provides insights into the microscopic examination of tissue for disease diagnosis, prognosis, and biomarker discovery. Particularly for cancer research, precise classification of histopathological images is the ultimate objective of the image analysis. Here, the time-frequency time-space long short-term memory network (TF-TS LSTM) developed for classification of time series is applied for classifying histopathological images. The deep learning is empowered by the use of sequential time-frequency and time-space features extracted from the images. Furthermore, unlike conventional classification practice, a strategy for class modeling is designed to leverage the learning power of the TF-TS LSTM. Tests on several datasets of histopathological images of haematoxylin-and-eosin and immunohistochemistry stains demonstrate the strong capability of the artificial intelligence (AI)-based approach for producing very accurate classification results. The proposed approach has the potential to be an AI tool for robust classification of histopathological images.

www.nature.com/scientificreports/ In this study, recurrent neural networks that learn time-frequency and time-space features for classification of time series or sequential data 15 is further developed for classifying histopathological images. The networks receive input as the combination of multiple features extracted in time-frequency and time-space domains of the time series transformed from the images. The class modeling is then designed for training the networks. Comparative results obtained from testing three public datasets of histopathological images, including haematoxylin-and-eosin (H&E) stained tissue images of colorectal cancer, H&E stained tissue images of heart failure, and immunohistochemistry (IHC) stained tissue images of rectal cancer, show the capability of achieving very accurate classification by the current approach.

Materials and methods
Image data. Three public databases are used in this study and descibed as follows.
H&E colorectal-cancer data. The colorectal-cancer (CRC) histology data used in this study were originally studied in 16 . The H&E stained tissue images of the CRC are publicly available at URL: http:// doi. org/ 10. 5281/ zenodo. 53169. The dataset consists of ten anonymized H&E stained CRC tissue slides. Tumors of both low grade and high grade were included in the dataset. The slides were first digitized, then contiguous tissue areas were manually annotated and tessellated to produce 625 non-overlapping tissue images of 150 × 150 pixels for each of 8 types of tissue, resulting in a set of 5000 images. There are 8 tissue types for classification, which are: (1) tumor epithelium (tumor), (2) simple stroma (stroma), (3) complex stroma (complex), (4) immune cells (lymphoid), (5) debris, (6) normal mucosal glands (mucosa), (7) adipose tissue (adipose), and (8) background (no tissue or empty). Figure 1 shows selected samples of the H&E stained CRC tissue images of the dataset.
H&E heart-failure data. The H&E stained heart-failure tissue image dataset, which includes left ventricular tissue from 209 patients was originally studied in 17 . The H&E stained tissue images of human heart failure are publicly available at the following URL: https:// idr. openm icros copy. org/ webcl ient. The dataset consists of two cohorts of patients: (1) heart failure (N = 94) and (2) without heart failure (N = 115). The heart-failure tissue was collected from patients with clinically diagnosed ischemic cardiomyopathy (N = 51) or idiopathic dilated cardiomyopathy (N = 43). The non-failing patients were organ donors without a history of heart failure. All tissue types were sectioned, stained, and scanned during the data acquisition. The whole slide image of each patient was down sampled to 5 × magnification, where eleven non-overlapping images considered as the regions of interest were extracted and the tissue border was manually refined. The total number of images for the heartfailure and non-heart-failure cohorts are 1034 and 1265, respectively, resulting in a dataset of 2299 images of 250 × 250 pixels. Figure 2 shows selected samples of H&E stained sections of heart-failure and non-heart-failure tissue types.  www.nature.com/scientificreports/ IHC rectal-cancer data. The data were obtained from a cohort of 143 patients with rectal cancer from the Swedish rectal cancer trial of preoperative radiotherapy (pRT) between 1987 and 1990 18 . All patients were performed locally curative resection. Among them, 77 patients received tumor resection alone, and 59 received pRT followed by surgical tumor resection. The pRT was given at a total dose of 25 Gy in 5 fractions over a median of 8 days (6-14 days) before the surgery. The surgical tumor resection was carried out in a median of 4 days (range 0-8 days) after the pRT. Samples were collected from biopsy ( n = 96 ), primary cancer or surgically resected tumor ( n = 136 ), adjacent normal mucosa ( n = 79 ), and distant normal mucosa ( n = 119 ). The distant normal mucosa was taken from proximal or distal margin (4-35 cm from the primary tumor) of the resected rectum and was histologically free from tumor. The adjacent normal mucosa sample was taken adjacent to the primary tumor and was histologically free from tumor. None of the patients received preoperative or adjuvant chemotherapy. The mean follow-up period or the patients was 107 months (range 0-309 months). The whole dataset has 235 images, where the size for each image is about 2500 × 2700 pixels. The image subsets consist of 40 and 14 images of biopsy without pRT having survival rate greater and less than 5 years, respectively; 32 and 11 images of biopsy with pRT having survival rate greater and less than 5 years, respectively; 54 and 25 images of tumor without pRT having survival rate greater and less than 5 years, respectively; and 36 and 23 images of tumor with pRT having survival rate greater and less than 5 years, respectively. Details about the procedure for the image extraction and the publicly available dataset for investigating the prediction and prognosis of DNp73 were described in 19 . Figure 3 shows selected samples of IHC stained images of the rectal cancer obtained for biopsy and primary tumor samples with and without pRT. Image vectorization. Because LSTM networks were designed for learning order dependence in time series or sequential data, the conversion of time-independent images into time series is necessary for extracting sequential features that will be used as input into the network. To extract features of the time series of histopathological images in time-frequency and time-space domains, the color (3D) images were first converted into grayscale (2D) images. Using the grayscale images, the vectorization of a 2D image or matrix is a linear transformation that converts the image into a column vector. Specifically, the vectorization of an M × N image I, denoted as J, is the MN × 1 column vector obtained by stacking the columns of the image I on top of one another, giving Thus, J results in a time series of image I, which is ready for feature extraction in time-frequency and time-space domains, which are described in the following sections.  www.nature.com/scientificreports/ or sequential data such as language modeling and speech-to-text translation 21,22 . The LSTM network was developed to solve the vanishing gradients problem encountered by the RNN. The essence of the LSTM design is the introduction of data-dependent configurations into the RNN cells to prevent the gradient of its objective function from vanishing during the data training 20 . As a result of the redesign, LSTM networks can be more robust and versatile than RNNs 23 .
As an extension of applications of the time-frequency and time-space LSTM (TF-TS LSTM) recently introduced in literature for classifying physiological signals 15 , images of histopathological tissue are vectorized into time series whose TF and TS features can be extracted for learning by the TF-TS LSTM. The extractions of the TF features that are the instantaneous frequency (IF) and spectral entropy (SE), and the TS features that are the fuzzy recurrence image entropy (FRIE) and fuzzy recurrence entropy (FRE) were described in detail in 15 . Details on how TF and TS features are used for learning by the TF-TS LSTM can also be found in 15 .
The procedure for implementing the TF-TS LSTM for image classification is as follows. The IF of a time series, which is the average of frequencies f over time instant t, is expressed as where P(t, f) is the power spectrum. Because Eq. (2) applies to infinitely long signal, the IF for a signal of a finite length needs to be numerically estimated. A method for estimating the IF and adopted herein was described in 15 .
The SE at time t is computed as where the probability at time t and frequency m is   www.nature.com/scientificreports/ in which f ∈ [0, fs/2] , and fs is the sampling frequency. Next, given a time series, an embedding dimension, and a time delay, the phase space of the corresponding dynamical system can be constructed and represented with a collection of vectors X = (x 1 , . . . , x N ) . Elements of an FRP are defined as 24 where µ(x i , x j ) ∈ [0, 1] is the fuzzy membership of similarity between x i and x j .
By partitioning X into c clusters using the fuzzy c-means algorithm 25 , an FRP has the following three properties: 1. Reflexivity: 2. Symmetry:

Transitivity:
The FRIE is then computed as where L is the number of gray levels of the FRP (represented by the fuzzy membership grades), and p l is the probability associated with gray level l.
Using the concept of the entropy of a fuzzy set 26 , the FRE of an FRP is calculated as 27

Class modeling and performance measures
Deep-learning based class modeling. In machine learning, methods for class modeling or one-class classification 28 aim to recognize samples of a particular class among all other class samples by focusing on the learning from a training set containing only the samples belonging to that class. This approach is different from conventional classification methods, which learn to differentiate two or more classes with training data containing samples from all the classes. In this study, the bi-LSTM network was used for the class modeling. The network learned samples from a particular class while creating an equal number of samples for other class(es) from a single sample for each of other class(es). This design was based on the capability of the bi-LSTM network to enhance its learning from multiple copies of the same sample. Not only this design of class modeling can enhance the learning capability of the network, it can also address the problem of data imbalance often encountered in the classification of histopathological images 29 . This deep-learning approach is referred to as image-LSTM.
The network layer was specified with an output size = 100, fully connected layer = number of classes, and, while the sigmoid function is used for the LSTM gate activation, the network ends with a fully connected layer, a softmax layer, and a classification output layer. Training options of the bi-LSTM were set as optimizer = ' Adam' (adaptive moment estimation), including L 2 regularization factor, maximum number of epochs = 80 for the 2-class classification of the H&E CRC, H&E heart-failure, and IHC rectal cancer data, and 180 for the 8-class classification of the H&E CRC data, minimum batch size = 150, initial learning rate = 0.01, and gradient threshold = 1.
To compute the instantaneous frequency and spectral entropy, the sampling frequency fs = 300 Hz, and frequency range = [0, fs/2]. In this study, the value for the sampling frequency is arbitrarily chosen to extract the time-frequency features of the transformed time series of the images, but commonly used for physiological signals (https:// physi onet. org/ chall enge/ 2017/). To compute the Shannon entropy and fuzzy entropy of an FRP, embedding dimension = 1, time delay = 1, and number of clusters = 3. For consistent input to the bi-LSTM, the length of the sequences of the time-space features (FRIE and FRE) was designed to match with that of the time-frequency (IF and SE) sequences by segmenting the transformed time series into the number of segments being equal to the length of the time-frequency sequences. Each segment of the time series was then sequentially processed to obtain the time-space sequences. Figure 4 shows the architecture and sequential procedure of the image-LSTM for classifying histopathological images.   For the prediction and prognosis of DNp73 in rectal cancer, survival rates of more and less than 5 years are defined as true positive rate (TPR) and true negative rate (TNR), respectively, which are calculated as

Results
To carry out the image data pre-processing for sequential feature extraction described earlier, because of the large size of the IHC slides of rectal cancer, these images were resized to 250 × 250 × 3 pixels to reduce the time for subsequent computations and network training. Figures 5, 6, and 7 show examples of the conversions of images into time series, time-frequency features, and time-space features using H&E CRC, H&E heart-failure, and IHC rectal cancer datasets, respectively. These features of the images were combined to constitute multi-dimensional sequences and used as the input into the bi-LSTM for learning and classification. Tables 1 and 2 show the classification results obtained from different methods using the H&E CRC and heartfailure data obtained as averages of ten runs of tenfold cross-validation (CV), H&E heart-failure data as averages of ten runs of twofold (for comparison with previously reported results) and threefold CVs (for comparison with previously reported results), and IHC rectal-cancer data as averages of ten runs of tenfold CV, respectively.
For the 2-classification (stroma and tumor) of the H&E CRC data, the image-LSTM provided perfect results in terms of accuracy, sensitivity, specificity, precision, and F 1 score, where applicable, and outperformed other  www.nature.com/scientificreports/ classification methods (support vector machines and histogram = SVM-histogram, and support vector machines and local binary patterns = SVM-LBP) using the same dataset reported in literature 16 .
Similarly, the image-LSTM outperformed the random forest (RF) and convolutional neural network (CNN) classifiers studied in 17 in testing the H&E heart-failure dataset using twofold and threefold cross-validation procedures for evaluating the performance of classifiers.
Once again, the image-LSTM provided excellent and better results (100%) than many pre-trained convolutional neural network (CNN) models recently reported 19 for classifying the biopsy and tumor samples with and without pRT in terms of accuracy, and survival rates > and ≤ 5 years.
To examine how quickly the image-LSTM accuracy was improving, and whether the network trainings could be overfitted with the training data, Fig. 8 shows the plots of the model training histories in terms of accuracy and loss using six different datasets for the tenfold CV. In general, if the loss keeps decreasing, the model could be overfitted; and when the loss becomes equal (decreasing to a point of stability), the model is considered either to be of a good fit or reaches a local minimum. By observing the plots for different datasets, the model quickly improved the accuracy and all the trainings did not result in overfitting.

Discussion
Regarding the classification of the H&E CRC data, there are much larger differences in the average accuracy obtained from the SVM-histogram (about 15% difference) and SVM-LPB (about 20% difference) classifiers for 2-class and 8-class classifications, where the latter case is lower than the former. The average accuracy provided by the image-LSTM classifier for the 8-class classification (99.96%) is almost the same as for the 2-class classification (100%). The better performance of the image-LSTM suggests that while the sequential features of the time series transformed from the original images helped increase the differentiation of the class properties, the formation of multi-dimensional sequential features allows the enhancement of the deep learning of long-term dependencies in more than one dimension. In fact, an important task of an LSTM is to decide what information to forget and keep from the cell state, and these two things are combined to make an update to the state. The time-frequency and time-space features enabled the updating process more effective, resulting in a better classifier.
The transformation of images into time series for extracting time-frequency and time-space features and class modeling using the bi-LSTM have shown the classification design effective as the results obtained from   16 80.80 ± n/m n/a n/a n/a n/a SVM-LBP 16 76.20 ± n/m n/a n/a n/a n/a image-LSTM 99.96 ± 0.08 n/a n/a n/a n/a H&E heart-failure: tenfold CV www.nature.com/scientificreports/ the image-LSTM outperformed the CNN, which is also a deep-learning model and used the direct input of the images, in both twofold and threefold cross-validations. H&E stained slides are used by pathologists for disease diagnosis and patient treatment. Images of cancer tissue are routinely examined by pathologists for cancer type identification and prognosis. The power of the precise classification of H&E stained images provided by the image-LSTM offers a useful AI tool to assist pathologists in carrying out timely analysis of large volumes of images.
The prediction and prognosis power of the protein marker DNp73 in rectal cancer can be discovered with the use of the image-LSTM for differentiating both biopsy and tumor samples of the rectal-cancer patients between with and without pRT. Recent results 19 showed the usefulness of several pre-trained CNN models for discovering the function of DNp73 in rectal cancer. In comparison with the pre-trained CNN models, the present classification results obtained from the image-LSTM not only confirm the recent findings but also provide a better AI tool for clinical decision making and accurate forecast of the future course of the cancer under pRT.
The high performance of the image-LSTM in classifying IHC slides can be very useful for timely biomarker discovery. Current practice in pathology relies on manual scoring of protein expression mainly based on colors to assess the capability of prediction and prognosis of the candidate protein as a biomarker. The classification accuracies based on the expressions of DNp73 with pRT and without pRT suggest the power of the present AI method with an implication that if the examination of a tissue taken from a rectal cancer patient is predicted to have a short survival ( ≤ 5 years), then the clinical decision would be to recommend the patient to be treated with pRT 30 .
While traditional IHC analysis did not provide any predictive and prognostic information of DNp73 expression in the rectal cancer patients without or with pRT 19 , both predictive and prognostic power of DNp73 expression can be discovered by the current AI approach. Such a discovery is very useful for optimal treatment and clinical decision making in pRT. The AI approach reported herein is not only useful for studying the protein expression in rectal cancer patients but can also be applied for discovering biomarkers in other types of cancer. A certain advantage discovered from the findings regarding the computer implementation of the AI is that accurate classification results can be achieved without the need for the manual extraction of regions of interest from whole-slide images. Such a freedom from the manual analysis is particularly useful for guaranteeing objective and reproducible results as well as alleviating time-consuming work to help life-science researchers focus on more important aspects in their study 30 .

Conclusion
An approach for classifying histopathological images by feeding the LSTM with sequential features extracted from time-frequency and time-space domains together with the class modeling has been presented and discussed in the foregoing sections. The key attribute is the capture of effective sequential features of texture-rich images in pathology that can leverage the power of the LSTM in learning sequential characteristics for pattern recognition. Classification results obtained from testing the new approach with both H&E and IHC image data www.nature.com/scientificreports/ outperformed other baseline methods. The image-LSTM appears to be a useful AI tool in the big data analysis of digital pathology for disease diagnosis, prognosis, and biomarker discovery, where the effective handling of big data can significantly contribute to modern healthcare systems 31 . The tool presented herein can be applied for classifying other types of microscope images of cells or tissues.

Data availability
The computer software written in MATLAB for the image-LSTM tool is publicly available at the author's personal homepage: https:// sites. google. com/ view/ tuan-d-pham/ codes.


extfieldsofstudy:
Physics
Computer Science

text:
Evaluation key generator of Multiple Asymmetric methods in Wireless Sensor Network (WSNs)

The wireless sensor network (WSNs) is constructed with a combination of the sensor nodes and the sink nodes. WSNs applied in many applications. The Security in WSNs is necessary because of limitations of storage, computing capacity, power use. Encryption is one of the most important tools utilized to supply the security services for WSNs. Asymmetric and symmetric encryption algorithm can be utilized in the WSNs structure to supply the security. The asymmetric key encryption algorithm gives a higher level of the security, but compared to symmetric key encryption, that causes extra sensor overhead. In this research KCMA method used to generate chains key of twelve experiments of the algorithms (ECC, RSA, ELGamal). These chains merged with hash function SHA2 and XOR. Diehard test was utilized in all experiments to evaluate the randomness of the secret key generated, and show it’s more security system. SHA2 was the best as compared with XOR. Also, the work evaluated the performance (time) for the system and throughput network.

Introduction
Sensor networks be composed of hundreds of self-organizing, low cost ,low power wireless nodes prevailed to observe the environment [1]. Sensor nodes capture information from an environment, process data and transfer it over radio signals. In sensitive WSNs applications like as enemy line monitoring or border area surveillance, the security protocols that supply confidential data transfer from sensor nodes to the base station should be utilized. However, the low processor and radio capabilities of the sensor nodes do not permit for the application of the traditional security protocols in the sensor networks [2]. The security is generally accomplished through encryption in WSNs. Encryption protects data through encryption and decryption. It composes of two classes: asymmetric and symmetric encryption, where one key is utilized to achieve the security process in the symmetric method and to various keys are utilized in the asymmetric method. Several encryption algorithms are suggested for security. Encryption algorithms work with key creation that guarantees message confidentiality.
The last years have known about the evolution of novel, more energy-effective encryption algorithms and give the same security threshold as traditional algorithms like as RSA [3,4]. The Elliptic Curve  [10] is one of these novel algorithms and is the most promising in terms of energy and time consuming, making it very attractive for encrypting data in WSNs. The ECC provides equal security with smaller key sizes, providing memory, computational and energy for restricted wireless devices [4] [9]. On the other hand, Elgamal was developed in 1984 by Taher ElGamal. It is an asymmetric key technique-based D-H key interchange. ELGAMAL encoding can be described on any periodic group G. Security depends on G problem related to the calculation of discrete logarithms [4] [9]. The rapid generalized encryption of long messages and the rate of data extending have the biggest advantages of this technique . The prime restriction in ELGAMAL is a random condition and a slower speed [5] [2].
The security is a very important issue in the WSN, especially with the development of science and computers ,the researchers attempted to provide security by generating chains secret key, the randomness of these keys generated are low, so in this work suggested Key Chain Multiple Asymmetric (KCMA) method to increase randomness and improve the confidentiality of the system. The key chains generate through twelve experiments evaluate with Diehard test. The time and throughput network for the chains generated was calculated in the Java programming language. Also, the complexity of the algorithms is illustrated.
The rest of the paper as follows: The WSNs are explained in the Section 1. The related work is highlighted in the Section 2. The proposed research is explained in the Section 3. The Measures Metrics and result is illustrated in the Section 4. The complexity of an algorithm is explained in the section 5. And, the Section 6 concludes the paper.

Related work:
The work related to researchers who work in this field: The search In [6], the research suggests an algorithm based on complementary subtraction 1 to represent numerical multiplication in numerical multiplication which provides a lower Hamming weight and significantly improves the computational efficiency of numerical multiplication .In this paper, used a hybrid cryptosystem asymmetric RSA and El-gamal algorithm , improve performance and the security will be relay on throughput and encryption time , decryption time [7]. This paper has shown increase the security system by using the protocol MKP which use two series of ECC algorithm to generate secret keys. It gives the authentication , since it has multiple key and work in parallel mode, that increase the throughput and the performance of the system ,and decrease the run-time, also it save the energy [8]. (Palanisamy and Chellappan 2011) suggested RSA -CRT technology is enhanced in secure routing and gives functional management of keys that use the best encryption / decryption to authenticate broadcast messages. The suggested technique will resist attacks on the WSN network layer and will minimize overhead for communications, and energy node consumption, and storage space [9].This paper [10] implemented ECC and RSA over Sun SPOTs . Its transmit the light reading, time and capacity of the available battery. Initially charged the battery, then computed time taken by every technique to discharge the battery to detect which technique exhausts minimal battery. In the same conditions RSA exhausts, power more than ECC. [11] in this paper used two system RSA and ECC,and discover that ECC is more suitable as compared to RSA, because of low CPU consuming, low memory use and shorter key size as compared to the RSA.

The Proposed Research:
The objective of the research is building a strong WSNs security system and improve the performance   To generate the elements in each chain, the asymmetric algorithms are used as in the following Eq.1 and Eq.2: Ai=PK1en (TK1, Ai-1) (1) 0< I <n Bi=PK2en (TK2, Bi-1) (2) 0< I <n Where: PK1en = the first algorithm, PK2en = the second algorithm. The keys are calculated using the hash function (f) as shown in the following Eq.3: (3) 1< i <n Where F represents hash function XOR or SHA2

Measured Metrics
This section demonstrated some important measures for the evaluation of the proposed KCMA method, this measure can be illustrated as follows:

Randomness test
The randomness of the output is one of the significant factors in measuring the security of any algorithm. The search implements diehard tests. Diehard (develop by Marsaglia) is a battery of statistical tests for measuring the quality of a random number generator.  [12] .These tests are p-values that dividing the range to three region : safe, doubt, and failure area. These regions can be defined by this following limits: 0.1 < p -value ≤ 0.25 or 0.75 ≤ p > < 0.9 lies in the Doubt region. 0 < p -value ≤0.1 or 0.9 ≤ p -value ≤ 1 located in the failure region . 0.25 < p -value < 0.75 located in a safe region [13].

Network throughput
Throughput is the amount of the data transferred in a unit of the time. Moreover, the node throughput is the rate of the data transported through the simulation time as described in Eq. (4). The throughput network is the throughput rate of all nodes on the number of the nodes as shown in Eq. (5) [14] T_N Data/Time… (4) Where T_N: the node's throughput, Data: the amount of the transferred data, and Time: the time simulation.

Randomness test
The test was performed on the twelve experiments of the asymmetric method described in the following table (1). The work file that size (11) MB was used as a measure to test the randomness of suggested research.
The randomness of the twelve experiments is illustrated in Figures below:

Time generating of the ECC methods
To calculate the time in five cases of the key (100,200,300,400,500), when using two ECC _XOR, it consumes the longer time of use it with SHA2, the time consumption of ECC_RSA and ECC_ELGamal when used them with XOR spent more time when used them with SHA2, but ECC_ELGamal is the least time consuming, as shown in Figure No. (16).

Time generating of the RSA methods
In the five cases of no. of the keys (100,200,300,400,500), we notice that the use of TwoRSA with XOR takes longer time than when used with SHA2, as well as ELGamal_RSA _ XOR consume time more than four ELGamal_RSA _ SHA2, as shown in the figure no. (17).

Throughput of the ECC based KCMA methods
From the figure (19), to measure throughput of five keys (100,200,300,400,500), the chart illustrates the high throughput network in the key (100), and network's throughput is less in the key number (500). The ECC_ELGAMAL_SHA2 is higher throughput, and TwoECC_XOR is a lower throughput network in this chart.    Figure 21. Throughput of the ELGamal based KCMA methods

Time Complexity
Time complexity is a mathematical complexity that describes the amount of the time it gets to run an algorithm. The complexity of the algorithm is classified according to the kind of function that appears in the large notation (O). For instance, a time complexity algorithm is linear time algorithm and a time complexity algorithm for constant is polynomial algorithm [15]. The problems of integer factorization of algorithm RSA permit algorithms which run in sub-exponential time [16]. according to [17]complexity of RSA refer to is (o(log(n)3)). With regard to El-Gamal encryption system is based on the discrete logarithm problem at Zp*,. It is necessary that different random integers k be utilized to encrypt various messages [18]. The time complexity of El-Gamal cryptosystem is (o(log(n)3) [17]. The complexity of ECC "difficult to break it" is equal to resolving problem of the discrete logarithm. Find discrete logarithm of one item in elliptic will not help locate the logarithm of any another item [19]. The time complexity of ECC is equal to(o√n) [20].

Conclusion
Security is one of the most difficult problems in WSNs. In this research, the KCMA method has been proposed to generate a key chain of asymmetrical algorithms. These chains generated from two different algorithms combined with hash function SHA2 and XOR, then it have been tested by Diehard Test, the results shown the randomness increased by using the SHA2, so the system is confidential and secure, as demonstrated by the test the highest security methods were (TwoRSA_SHA2, TwoELGamal_SHA2, RSA_ELGAMAL_SHA2, ECC_RSA_SHA2, ECC_ELGamal_sha2). The time was calculated in the Java language for the generated keys and, shown that the method which used less time consuming were (ECC_RSA_SHA2, ECC_ELGamal_ SHA2, ELGamal_RSA_SHA2, ECC_RSA_SHA2, TwoELGamal_XOR). SHA2 is the best in terms of security and time consumption. The complexity of the RSA and Elgamal were o(log(n)3, and ECC is o√n.


A Classification Method of Acupoints and Non-acupoints based on Traditional Features and Wavelet Features

Meridians and acupoints are the basis of TCM theory and play an important role in disease diagnosis and acupuncture treatment. There are still many problems in current research on electrical signals of acupoints. On the one hand, most of the studies did not consider the integrity of the meridian, but only based on a few acupoints. On the other hand, the lack of targeted feature extraction and classification methods leads to unsatisfactory classification results. Considering the above problems, a method combining traditional features and wavelet features is proposed to classify acupoints and non-acupoints. Based on the integrity of the meridians, we first collect the body surface electrical signals of some acupoints and non-acupoints on the twelve meridians of the human body, and then extract traditional and wavelet features from the measured signals. Finally, SVM and XGBoost are used to classify acupoints and non-acupoints respectively. The experimental results show that this method can effectively improve the classification performance of acupoints and non-acupoints, and for the feature vectors constructed in this paper, XGBoost has better classification capabilities.

Introduction
In recent years, the problems existing in the western medicine therapy, such as chemical drug resistance and increase in the cost of medicare, cannot be neglected any more. Traditional Chinese Medicine (TMC) has attracted more and more attention because of its safety, convenience and effectiveness. However, compared with the western medicine, the application of TCM is constrained in many areas due to the lack of objective and scientific basis in meridian theory [1] and the absence of a clear and independent anatomical structure.
Acupuncture points are the basis of traditional Chinese medicine theory, and their states can reflect the physical health of the human body. Therefore, conducting research from the perspective of electricity, analyzing their electrical characteristics, and distinguishing them from non-acupuncture points, can help to achieve precise positioning of acupuncture points. And it can also provide similar research ideas for disease diagnosis and treatment through acupoint electrical signals, which has certain significance for supporting the theory of Chinese medicine.
The current research on the electrical characteristics of meridian points is mainly divided into two directions: one is to find out the different characteristics between acupoints and non-acupoints under the condition of applying electrical stimulation, and some results have been achieved. For example, a large number of experimental results show that the acupuncture points have low impedance and high capacitance characteristics [2]. The other is to study the potential characteristics of acupoints without stimulus, and it has been concluded that acupuncture points have higher potential characteristics than non-acupuncture points [3]. Compared with the method of applying electrical stimulation,The study of acupoint potential without stimulation directly collects the electrical signals spontaneously generated by the human body for analysis, which avoids the influence of external stimulation, so it may be closer to the essence of the electrical characteristics of acupoints [4]. However, there are still many problems in current research on acupoint potential, such as: 1) Most studies are only based on the potential characteristics of a few acupoints and nonacupuncture points on one meridian. They have not studied the similarities and differences between acupuncture points and non-acupuncture points on multiple meridians based on the integrity of the meridian.
2) The lack of targeted feature extraction and classification methods results in low classification accuracy.
3) There are relatively few studies on the electrical signals of acupuncture points and nonacupuncture points without stimulation.
In order to solve the above problems, this paper collects the electrical signals of some acupoints and non-acupuncture points on the twelve meridians of the human body, based on the integrity of the meridian, considers both the global and local features of the signals, and finally uses machine learning methods to classify acupoints and non-acdupoints.
The rest of this article is organized as follows. Section 2 briefly introduces the overall plan of this research. Section 3 mainly introduces the data processing process. Section 4 shows the experiments and results. Conclusion and future work are given in Section 5.

The overall plan
The overall plan of this article is shown in Figure 1. Firstly, we collect the raw signals of acupoints and non-acupoints on the twelve meridians of the human body, and then denoise the raw signals. Next, we extract traditional features and wavelet features from the de-noised signals. Traditional features consider the overall information of the signals and describe the overall characteristics of the signals, while wavelet features consider the local information of the signals and describe the local characteristics of different sub-band signals. Finally, machine learning methods are used to classify acupoints and nonacupoints.  3 points and some non-meridian non-acupuncture points 2cm near the acupuncture points on the twelve meridians were collected sequentially. All measured points are on the left body of the subject. At each measured point, we collected a 1000hz,150s-long signal . A group of subjects includes 7 healthy volunteers (20-29 years old, 3 males, 4 females). For a single subject, we obtained the raw signals of 47 acupuncture points and 38 non-acupuncture points.

Data pre-processing
For a single 150s long signal, we cut it every 3s as an acupoint or non-acupoint signal data. Among the 7 subjects, we selected the collected signals of 6 subjects (3 males, 3 females) as the training set data, and took the collected signals of the last subject as the test set data. In the training set data and the test set data, the total amount of acupuncture point and non-acupoint point samples is balanced.
We first observed the frequency domain of the collected signal and found that the 50hz power frequency and its 100hz harmonic interference in the collected signal frequency domain were relatively large, so we used the notch filter to remove the 50hz and 100hz power frequency interference [5].For the signal whose power frequency interference has been removed, we used wavelet decomposition filtering techniques [6] to remove other noises. As show in figure 2, comparing with the original signal, the denoised signal eliminate useless noise on the basis of retaining the useful signal to the maximum extent, which is very beneficial for feature extraction in the following work.

Traditional features
Traditional features include time domain features, frequency domain features, and nonlinear features. For a denoising signal, this paper extracts 42-dimensional traditional features.

Time domain features.
We perform various operations on the amplitude of the signals to extract the time domain features which can intuitively measure the change of the measured signal's amplitude and shape .The statistical characteristics used in this paper are: count value, average value, maximum, minimum, range, variance, skewness, kurtosis, etc. The method for calculating the variance is given by equation (1), where X i is the value at each time, N is the length of the sequence and 〈x〉 is the mean of the sequence.
The calculation method of the kurtosis shown in equation (2) can characterize the sharpness of the data distribution.

Frequency domain features.
Due to the randomness of the acupoint signal, the time-domain analysis has a certain instability, so we further obtain the signal's frequency-domain characteristics, which further reflect the nature of the body surface electrical signal.
In the frequency domain analysis, the measured data need to convert from the time domain to the frequency domain by FFT(Fast Fourier Transform), and then extract the corresponding features. The features used are: Shannon Entropy, Mean power, Total power, signal-to-noise ratio(SNR), etc. The Shannon entropy can represent the degree of uncertainty of the sequence. The calculation method of total power shown in equation (3) can represent the sum of the power spectral density,where p i is the power spectral density.
(3) The calculation method of signal-to-noise ratio shown in equation (4) can be used to characterize the proportion of valid information in the signals, where s is the signal's total power and n is the noise's total power. SNR = s/n (4)

Nonlinear domain features.
These features are extracted based on the nonlinear analytical method. The features in this paper include coefficient of variation(COV), median absolute deviation(MAD), autocorrelation coefficient(AC), the number of turning points (NT) [7], etc. The calculation method of COV shown in equation (5) can characterize the degree of dispersion of the signal, where std(X) is the variance of the sequence. COV = std(X)/〈x〉 (5) NT counts the number of changes in the sign of the slope, in other words, the number of signal peaks, and the calculation method is shown in equation (9), where u(x) indicates a unit-step function.

Wavelet features
Wavelet analysis transforms a one-dimensional signal in the time domain into a two-dimensional timescale space. The time-domain signal is projected and decomposed at multiple scales in the wavelet transform domain to obtain signals of different sub-bands. Different sub-band signals contain different local information and describe the dimensions of the signal more finely [8]. The basic principles of discrete wavelet transform are as follows: Assuming that there is a continuous signal f(t), we perform discrete wavelet transform on it. The definitions of discrete wavelet transform and inverse transform are shown in equation (7) and equation (8): where ψ j,k (t)=2 -j/2 ψ(2 -j/2 t-k) is a wavelet sequence, and ψ(t) is a wavelet basis function that satisfies certain conditions. f j (t) represents the component of the signal f(t) at a certain scale (2 j ). j and k respectively represent frequency resolution and time offset.
By using the mallat algorithm, we decompose the signal f(t) in a finite layer, then f(t) can be expressed as equation (9): (9) In equation (9), L is the number of decomposition levels, D j is the detail component at different scales, and A L is the low-pass approximation component,more specific equation, A L =f L A (n) and D j =f j D (n). The entire frequency band of the signal is divided into multiple sub-bands. The approximation coefficients and wavelet layer coefficients corresponding to these sub-bands are cA L ,cD L ,cD L-1 ,…,cD 1 .Let us define f(n) as the sampled signal of continuous signal f(t), then f(n) can be regarded as an approximate value when the scale j=1, that is, A 0 (n)=f(n). As shown in Figure   AIACT 5 3, the discrete signal is decomposed by wavelet with scale j=2,3,4, and we get the approximation coefficient cA 3 and wavelet coefficients cD 3 ,cD 2 ,cD 1 . Due to the lack of prior knowledge of acupoint signals, a thorough decomposition of acupoint signals may be more helpful to extract useful local information. In this paper, we use the Daubechies (DB7) wavelet to decompose the signals in 7 levels (j=8), and get the wavelet decomposition coefficient vector C. The structure of the coefficient vector C is as follows: C = {cA 7 ,cD 7 ,cD 6 ,cD 5 ,cD 4 ,cD 3 ,cD 2 ,cD 1 } where cA 7 is the approximation coefficient, and cD 7 ,cD 6 ,…,cD 1 are the wavelet coefficients.
For the approximation coefficient and the wavelet coefficients of each band, six statistical features which include maximum, minimum, mean, median, standard deviation, and average energy are extracted respectively. Therefore, an original signal can obtain 48-dimensional wavelet features [9]. These wavelet features describe different local information of the original signals.

Experiments and Results
In this paper, we obtained 7112 1000hz 3s-long signals from the original signals of 6 subjects, which include 3920 acupoint signals and 3192 non-meridian non-acupoint signals. These signals are used as training set data.
In order to evaluate the classification ability of the model more accurately. We randomly selected 12 acupuncture points and 12 non-meridian non-acupuncture point signals from all the measured signals of the seventh subject as a test data set. We randomly selected 11 times and got 11 different test data sets. Each test data set has 672 1000hz 3s-long original signals ,which include 336 acupoint signals and 336 non-meridian non-acupoint signals.
We take traditional features as the first type of features, and the combination of traditional features and wavelet features as the second type of features.
Firstly, we extract the first type of features from training data and test data. In the training stage, we use the training feature data to train SVM based on the RBF kernel [10] and XGBoost based on the CART regression tree [11], optimize the parameters by the grid search method ,and finally save the trained model. After that, we use 11 test data sets that have the first type of features to verify the classification effects of different models. F1-score is used as the evaluation standard for the classification ability of our models. Then we extract the second type of features from the training data and test data, and repeat the above process. As we have seen in Figure 4 and Table 1, from the perspective of the features, when using the second type of features that combine traditional features and wavelet features, the average F1-score of SVM in the 11 test sets is 0.7373, which has an increase of 0.0392 compared to using only traditional features, the average F1-score of XGBoost in the 11 test sets is 0.7728, which has an increase of 0.0292 compared to only using traditional features. The experimental results show that, compared with the first type of feature vectors that only extract the global features of the signals, the second type of feature vectors integrate the global and local features of the signals, have higher separability, and improve the model's ability to classify acupoint and non-acupoint potential signals. From the perspective of the models, using the feature vectors constructed in this paper, the average F1-score of the SVM and XGBoost classification models on the 11 test sets are all above 0.73, which can achieve the purpose of classification for acupoints and non-acupoints well. And when using the same type of feature vectors proposed in this paper, the classification performance of the XGBoost is better than SVM, so it has a better generalization ability.

Feature importance
In order to further explore the contribution of traditional features and wavelet features in the second type of features, We calculate the importance of all the features in the XGBoost model , shown in Table2, which show Top10,Top20,Top40,Top60 important features proportion. We can find that the traditional features are effective and they represent the global features of the signals. Meanwhile, the wavelet features are also indispensable, they represent the local characteristics of the signals. Combining global features with local features can describe the regular differences between the feature spaces of acupoints and non-acupoint signals more comprehensively, thereby