Throughout the history of mankind, from the ancient Egyptian and ancient Spartan scytale to the Caesar cypher and the Enigma machine in World War II, the development of cryptography guarantees the security of information storage and transmission. Those kinds of cryptography belong to symmetric-key cryptography, which utilizes the same cryptographic keys for encryption and decryption. The weakness of those strategies is that once the key tables are exposed, the cipher failed immediately. Apart from the symmetric-key cryptography, an innovation of asymmetric cryptography was proposed in the 1970s. Rivest-Shamir-Adleman (RSA) is one of the famous and successful examples of the asymmetric cryptography, which is widely used in information security¹. The security of the RSA cryptosystem relies on the difficulty of taking discrete logarithms for finding large prime factors of a given product. Although the RSA algorithm provides significant asymmetry between the intended and unintended uses of encoded information, the random number-based key generation for this algorithm is still potentially vulnerable (e.g., the forward secrecy issue). Threatened by the multifarious brute-force attack, cyber-attack, or Trojan implantation from unauthenticated parties, the random numbers in the security system should be intrinsically difficult to replace and predict^2⇓-4. Thus, the random number taken in the cryptosystem should be truly random or chaotic in nature. From the above considerations, implementing the physical or hardware primitives into cryptosystem could be a new attempt to enhance the robustness and complexity to guarantee the information security in the current big data era^5⇓⇓-8.

Memristor, a key component towards in-memory logic and neuromorphic computing, has realized versatile functions based on its cross-pointed structure⁹. However, for practical applications, the memristor's intrinsic stochasticity on its switching voltages, times and energies, as well as the fluctuations of the resistance state over time should be carefully considered¹⁰. Although many efforts have been made to control the CF location and growth, e.g. nanocone shape Ag electrode¹¹ or graphene as insertion layer for blocking excessive atoms/ions injection into the dielectric,^{12⇓⇓⇓-16} the intrinsic memristive stochasticity is in turn a naturally valueble entropy source for cryptographic research and applications^{17⇓⇓⇓⇓-22}. In the past few years, especially very recently, the use of the intrinsic memristive variability has aroused world-wide scientific attention in true random number generator (TRNG) as well as cryptographic key generation. In 2015, Balatti et al. constructed a TRNG using the statistical variability of the Set voltage in an oxide-based resistive switching²³. In 2016, Wei et al. at Panasonic developed TRNG based on memristive current fluctuation²⁴. In 2017, Jiang et al. utilized the volatile switching in a Ag:SiO₂ diffusive memristor to generate high-speed random bits²⁵. Following the work of the Jiang group, Woo et al. constructed a TRNG based on Pt/HfO₂/TiN memristor in 2018²⁶ and implemented a nonlinear-feedback shift register into the same system in 2020²⁷. In 2021, the same Woo group reported high-speed TRNG based on Cu_xTe_1-x diffusive memristor²⁸. And in the same year, Kim et al. proposed a NbO_x Mott memristor-based TRNG based on its intrinsic stochastic oscillation²⁹. In 2021, Wen et al. realized advanced data encryption using h-BN based memristor based on its intrinsic random telegraph noise (RTN)³⁰. In 2022, Liu et al. demonstrated the bi-modes entropy sources in a Bi₂O₂Se-based memristive TRNG, including statistical variability of Set voltages and current fluctuation over time³¹. Although many successful TRNGs based on the intrinsic memristive stochasticity have been reported, comprehensive studies of the intrinsic memristive stochasticity are still yet to be carried out.

In this study, the memristive spatiotemporal variability within graphene CBRAM has been synergistically investigated based on kinetic Monte Carlo (KMC) simulation, multiple Artificial Intelligence techniques, kernel density map and time constant analysis, which verified the inherent randomness of the memristive stochasticity. The random numbers based on hardware primitives are also passed through the Hamming distance calculation to verify their randomness, uniqueness and reproducibility. For the device architecture, graphene-based conductive bridging RAM (CBRAM) was chosen in the current study due to the dimensional anisotropic characteristics of graphene electrodes: high electrical and thermal conductivity in the plane direction, and low electrical and thermal conductivity in the vertical direction³². This characteristics could facilitate the transportation between device and peripheral circuit as well as inhibit the power consumption during the resistive switching. To analyze the spatiotemporal variability, three types of artificial intelligence (AI) technologies were employed including principal component analysis, long short-term memory and Shapley value. Comparing with the previously reported in-situ material detection or statical calculation, the strategy employed in the current study based on AI technologies provides the D2D correlation and classification analysis and the C2C correlation in long time series. Moreover, integrating the random numbers based on memristive variation into a Rivest-Shamir-Adleman (RSA) cryptosystem has been demonstrated in detail.

Unlike the metal-oxide-semiconductor field-effect transistors (MOSFET), whose working mechanism relies on a large number of electrons and holes, memristor, in contrast, realizes its resistive switching and information storage only via a limited number of ions within the switching medium³³. Fabrication inhomogeneity (e.g. the poly-grain nature of the graphene electrodes, the polymer-assisted transfer/lithography processes, and the subsequently nonuniformity deposition), and the initial formation of conduction filaments as well as its intrinsic stochastic switching mechanism, lead to the memristive device-to-device variability (DDV). Note that the device structure in the current study can be seen in the Fig. 1a and the fabrication details can be found in the Experiments section and Fig. S2 in the Supplementary Materials. A typical resistive switching curve of graphene-based CBRAM have been shown in Fig. 1b. For the Set process, the active electrode (silver in this study) gets oxidized into Ag⁺ and dissolved into the electrolyte, which is then reduced at the inert electrode surface and piled up as a metallic filament to bridge the active electrode during the Set bias voltage, during which the device is switched to the low resistance state (LRS). When an opposite polarity Reset voltage is applied, the conductive filament starts to be ruptured and the device is switched back to the high resistance state (HRS), which is the Reset process. Considering the differences of fabrications, initial filament formation and variation of repeated cycles and current fluctuation over time, the filament revolution (from Ag to Ag+ and from Ag+ to Ag) can be in highly stochastic, as shows in Fig. 1c. Since the silver cations transport more easily in a defect-rich region with lowest electrochemical potentials, the variability of the CF shape will be formed due to the inhomogeneous nature of the AlO_x³⁴. Briefly, the AlO_x amorphous electrolyte plays at least 3 roles for the variability during forming process and subsequent multiple cycles^35⇓-37: 1. noise source for current fluctuation: the electrolyte may include redox reactions of the Ag atoms or ions at the filament surface, or resulting atomic reversible rearrangements at the metastable positions or the defects sites, which may cause current fluctuation over time; 2. solubility: the CF and sub-filaments (if existed) will be dissolved at various degrees during the gradual Reset operation, which determined the HRS in the following cycle; 3. thermal conductivity: during repeated bias sweeping back and forth, a large amount of Joule heat is generated and accumulated within the electrolyte; Ag nanoparticles will merge into larger clusters to minimize the interfacial energy. As a result, not only the shape of the CF, but also the device performance will be influenced, because both the drift and diffusion current of the Set and Reset processes are thermally accelerated; in other words, whether the heat can be promptly dissipated with the electrolyte is also a stochastic factor to the device switching characterizations. Moreover, considering the weak thermal conductivity of the van der Waals interface of graphene electrode, the memristive variation will be further complicated. To analysis the CCV statistically, conventionally, a large number of IV curves could be analysis via Weibull distribution, time series analysis (TSA) or multivariate strategy to forecast the Set and Reset jointly³⁸. However, they neglect the physical correlation between the cycles. To consider the correlation between cycles, an iterative recurrent neural network based Long Short-Term Memory (LSTM) has been utilized to analyse the CCV of oxide-based memristor successfully. In this work, we extracted the physical parameters from memristive measurements and utilized principal component analysis (PCA) to analyse the DDV and LSTM to analyse the CCV. Moreover, the TTV is also analysed based on physical models, where the TTV under different voltage biases are shown in Fig. 1d1-d3. The device array and the device architecture have been shown in Fig. 1e-g via scanning electron microscope (SEM) observation. To inspect the device vertical structure and the thickness of the electrolyte layer, the transmission electron microscopy (TEM) image is shown in Fig. 1f. Moreover, elemental distribution of the TEM image with energy dispersive X-ray spectroscopy (EDX) including Ag, Al, C, O and Si has been exhibited in Fig. S3a-g in Supplementary Materials. Raman spectrum detecting has been also carried out to verify the high crystallinity and low defect density of the graphene electrode, as shown in Fig. S3h in Supplementary Materials.

To provide more physical understanding and better visualisation of the DDV and CCV during the device operations, a kinetic Monte Carlo (KMC) Simulation is carried out, and described in this section^39,40 (more calculation details could be found in Experiment section)⁴¹. As shown in Fig. 2a, a series of physical processes are considered in the simulation including anode oxidation (p1), cation hopping within the electrolyte (p2), cation reduction and deposition on cathode (p3), electron emission (p4), electron hopping within electrolyte (p5), electron adsorption (p6), cation and electron reduction (p7). The Fig. 2b illustrate the migration energy barrier of cation and electron, where γ and δ are factors of barrier lowing and a is the effective hopping distance of particles⁴². For reduction process, kq2/a is the barrier lowing effect due to Coulombic attraction. As shown in Fig. 2c, the Forming process of 4 devices are under same operation and the same physical structure, the filament morphology is distinct from each other (4 more Forming simulation results can be found in Supplementary Materials, Fig. S1a). The Forming process determined the initial conductive filament distribution and locations, as well as the distinct DDV. Considering the poly-crystal nature of graphene and the inevitable polymer residues during the transferring process, the discrepancy of graphene lattice and surface could also improve the complexity the DDV. After the Forming process, a sequential resistive switching has been carried out, as shown in Fig. 2d. Obviously, around the location of the initial filament, there is a higher possibility to form new filaments. Subsequently, during cycle-to-cycle Set and Reset, the filament forms and ruptures in different locations and morphologies. Extended with the Fig. 2d, 8 more simulation results could be found in Supplementary Materials, Fig. S1b. The simulation results are consistent with the recent report on in situ observation of a planar Pt/PEO-Li⁺/Ag memristor devices⁴³. Based on the direct observation of the CF growth revolution under different devices and repeated cycles, it is evident that the shape, formation position as well as the sub-filaments are in highly variable. The amorphous nature of the electrolyte will result in device to device variation in switching performances. And the preferential formation of CF does not prevent reformation of filaments in new locations, which reveals a high variability from cycle to cycle. As an alternative to the time-consuming in-situ characterization techniques, KMC simulation could provide more results based on the physical understandings and software calculations⁴⁴. However, moving towards the application of cryptographic key generation, further study to prove that the memristive variation is really random is necessary.

Traditional analysis tool presumes the independency between the physical parameters. Moreover, only few parameters, such as Set and Reset voltages are concerned. But based on the physical understanding of the memristor, the CCV is connected from cycle to cycle, e.g., the tunnelling gap formed in the previous filament rupture will affect the new filament formation. In turn, the shape, size, component and element gradients in the new formed filament will also influence the filament rupture in the subsequently cycles⁴⁵. Scaling to a broader spatiotemporal range, the DDV not only includes the fabrication differences but also the initial formation of the conduction filaments, as well as CCV. TTV is also based on the resistance state, which is also affected by CCV (in other words, the CCV determined the TTV). Thus, utilizing the spatiotemporal variability as hardware primitive, a critical question is that, are DDV and CCV independent of each other or are they physically correlated? Once the attacker obtains their physical correlations, the cryptosystem based on the memristive primitives could be easily encoded.

Although the KMC visualised the DDV and CCV results based on simulation, it cannot tell whether those results are correlated or not. To further analyse the DDV and CCV, artificial-intelligence-based techniques could be a promising strategy⁴⁶. Eight physical parameters have been extracted from IV curves including the Set and Reset operations (Set Voltage, Current before Set, Reset Voltage and Current, Current and Voltage after Reset) as well as the high resistance state (HRS) and low resistance state (LRS). The HRS and LRS are extracted at the voltage of ± 0.1 V, where the 0.1 V is for the HRS and the -0.1 V is for the LRS. The current and voltage of Set and Reset are extracted from the points after Set and Reset behaviours. The current and voltage after Set and Reset are extracted from the points after Set and Reset. Note that the Set processes are instantaneous, which is easy for extraction. And the Reset processes are gradual, where the lowest current is defined as the switching point for data collection. Before being fed into the following machine learning techniques, all current parameters were subjected to logarithm calculation for feature scaling. To analyse the DDV, we utilized Principal Component Analysis (PCA) to deduce whether the physical parameters from the identical device could be classified into an identical group⁴⁷. The PCA is to project a high-dimensional node into a low-dimensional space, and the low-dimensional space retains most of the properties of the high-dimensional space. By transforming a sample with n feature spaces into a sample with k feature spaces, where k must be less than n. Such transformation only allows linear transformations. In machine learning, PCA is a frequently used method of feature extraction for dimension reduction, of which the purpose is to reduce the number of dimensions of input data to improve the learning performance, but the overall prediction performance will not be degraded and even become better. 8 devices, each with 5 successive cycles have been taken into calculation. As shown in Fig. 3a, no clear classification can be obtained, which could be a solid evidence for the independency between DDV.

After the DDV analysis, CCV has been analysed with LSTM technique. The LSTM cell is shown in Fig. 3b. Based on its input gate, forget gate, activation gate, and output gate, the LSTM cell can update the cell state, remove ignorable information and propagate useful information to the next cell via feedback optimization and iterative calculation⁴⁸. LSTM is a type of recurrent neural network that allows the outputs at a previous time step (physical parameters in the previous cycle in the memristive calculation) to be conjunctively used as co-inputs at the following time step⁴⁹. Thus, the cycle-to-cycle correlation can be considered in the pointwise multiplication. In this scenario, to predict the Set voltage in a cycle (whether larger or smaller than the threshold voltage), such a recurrent neural network can process not only single data points (the rest of 7 physical parameters), but also entire sequences of 256 cycles. The feature extraction method is shown in Fig. 3b2 and the extracted physical features of those 256 cycles have been shown in Fig. 3c-j. Those extracted physical features are in highly random distribution, which is consistent with the KMC simulation results. To visualise the contribution of those physical parameters to the Set voltages, Shapley value plots were carried out to exhibit the correlation and contribution of CCV. The Shapley value plot is originated from game theory and local explanations, which could concisely leverage and visualise all the features’ contributions into the model's prediction. As shown in Fig. 3k-l, the “Reset Current” and “voltage after Reset” are located at the highest value of the latitude variables, indicating their biggest contribution to the prediction of the Set voltages. There two features are related to the start and the end of the Reset transition process. Although the “Set voltages” belong to the Set process, the physical features in the Reset process provide much bigger contribution and variation, which is due to its higher randomness in nature. Considering the heat accumulation effects, both the drift and diffusion processes will be influenced⁵⁰. Thus, the gradual Reset process is more complex in nature and also contributes the highest impact on the Set voltages. The “current after Reset” is also an important feature of Reset, but with a much lower contribution to the Set. This result could be attributed to its lower variability compared with the voltage. In other words, the variation of the current amplitude is more dominant during the Reset process. Subsequently, the HRS and LRS also play important roles in “Set voltages”. The HRS is related to the filament gap formed during the previous rupture and the LRS is related to the size, shape, conductivity, and ion gradient of the current filament. These calculation results are different from our recent study, a LSTM analysis of the CCV in a Bi₂O₂Se-based memristor, which may be due to the different switching nature between redox-based and conductive-bridge-based memristors^31,40. The SHAP values exhibiting not only the contribution but also the correlation between input value to the output are visualised. The red dots and blue dots located at the right x-coordinate line indicated a positive/negative contribution to the output value (the Set voltage), and vice versa. As shown in the Fig. 3 a, the red and the blue dots are synergistically overlapped for all the parameters. This result elucidated the inherent stochasticity of the CCV (because the 7 physical features synergistically influence the Set voltage, which is hard to be predicted physically), machine learning attack resilient and physically unpredictable features as hardware entropy source for the random crypto key generation⁵¹. The PCA techniques could visualise the classification of physical parameters of DDV and the LSTM with SHAP techniques could analyse in-depth the CCV within single device. Combining these AI techniques could clarify the source of the electrical variation, both for DDV or CCV, which could be a general analysis solution for other device platforms.

Based on the machine learning analysis, the DDV and CCV could be served as entropy sources for cryptography key generation. The random telegraph noise (RTN) is a stochastic current fluctuation between two resistance states. A 0.1 V constant voltage stress (CVS) in 100 Hz with different back gate voltage including 0 V, -0.1 V and -0.2 V, was applied on the graphene based memristor on its HRS state, due to its higher fluctuation than the LRS. Fig. 4a shows the power spectral density (PSD) of three kinds of RTN (0 V_BG, -0.1 V_BG and -0.2 V_BG), evidencing a 1/f behaviour. In a recent study of the RTN based on scanning tunnelling microscope point contact devices, the atomic fluctuation is verified as the dominant noise source⁵². The noise solely originated from dynamical defects, where atoms are fluctuating between metastable positions and driven by temperature-activated Langevin dynamics⁵³. MATLAB script was carried out to derive the flipping time of the three kinds of RTN signals. Most of the rising and falling transition time located at the region of 0.01 s symmetrically, indicating the same flipping time of rising and falling, as shown in Fig. 4b1 and b2. The time lag plot (TLP) of the three kinds of RTN signals in Fig. 4c-e, shows the distinctive two resistance states, where the left bottom corner is the lower resistance state δ₀ and the top right corner is the higher resistance state δ₁. Note that there is no transition state in the current study (no points located in the right bottom and top left corners)⁵⁴. Based on the symmetry distribution of the two states and equal transition time between the two states, the RTN signals provide equal probability for upper and lower current levels. In this scenario, the RTN signals (as shown in Fig. 4f1) could be utilized into random number generation via a binarized algorithm as shown in Fig. 4f2. The Fig. 4h-j exhibited the random number generation from the three kinds of RTN signals. Fig. 4g shows the D flip-flop and comparator-based peripheral circuit to extract and store the RTN signals into the binarized type. More RTN results under 20 Hz, 10 Hz and 5 Hz with 0 V, -0.1 V and -0.2 V back gate voltages have been shown in Fig. S4a-i. The lower frequency-based measurements do not exhibit adequate random bi-state bits, which may be due to their lower detecting rate comparing with the flipping rate of the RTN.⁵⁵

To generate cryptographic key, the intrinsic stochasticity of the DDV, CCV and TTV have been analysed on the physical level. Before harnessing those hardware-based primitives into the cryptography algorithms, Hamming Weight, inter-Hamming Distance (inter-HD) and intra-Hamming Distance (intra-HD) have been calculated to verify their randomness⁵⁶. The Hamming weight is the number of non-zero symbol of a string, also equivalent of the Hamming distance with an all-zero string in the same length. As shown in Fig. 5a, the Hamming Weight distribution could be well fitted with Gaussian where the mean value is approximately in 50%, indicating an equal probability for the generation of “0” and “1”. The inter-HD is defined as the uniqueness between different devices and the intra-HD described the intrinsic variation from the same device. The inter-HD and intra-HD are all around 50% as shown in Fig. 5b-c, indicating the randomness of the bit generation (the calculation details could be found in experimental section). Note that the intra-HD is derived from 528 cycles from same device, where each key has 24 bits; and the inter-HD is derived from 26 devices and each device contains 48 cycles (in other words, 26 keys from different devices, each key has 48 bits). Combining the DDV, CCV and TTV could enhance the robustness of the cryptography keys, which include spatiotemporal variability together, as shown in Fig. 5d⁵⁷. Those key lengths could be enough for the lightweight cryptography in embedded system or sensor networks. Note that the combination with periphery circuit could further improve the performance of the TRNG, e.g., the two-memristor-based structure could guarantee random number generation when one of the memristor are at fault stuck²⁶; a combination of nonlinear-feedback shift register (NLFSR) could further improve the complexity and generation rate of the TRNG²⁷. Fig. S6 exhibits three types of NLFSR, including Geffe generator, Massey-Rueppel multispeed generator, and Beth-Piper stop-and-go generator.

Implementation of the hardware primitives into cryptosystem can be a new fashion in both cryptology and electronic area. Here an example of hardware random numbers based asymmetric cryptosystem is shown in Fig. 5e. And the crypto key generation could be based on the RSA algorithm as shown in Fig. 5f. Note that, as comparison, a symmetric cryptosystem including text encryption and image encryption is shown in Fig. S5a1-c4 and Note II. Through a bit XOR method⁵⁸, a message of “OLYMPIC” could be encrypted as “,/;.3* ” and an image of the device structure of current study could be successfully encrypted (if the RGB histogram can provide useful information, attackers could derive it and obtain information of the original image). Back to the RSA cryptosystem, through the prime number selector (as shown in Table S1 in Supplementary Materials), the random numbers have been transformed into p and q values, which are the key random numbers in the RSA system. As a result, the public key (N, E) and the private key (N, D) could be derived based on the hardware primitives. Then the demonstration of the private key and public key in information security is shown in Fig. 5e. Supposing that Bob wants to say “Hi Alice” to Alice, there are two paths to guaranteeing the information security with the private keys and public keys. In path 1, Bob encrypted the message by using the Alice's public key and sent to Alice; then Alice utilized her own private key to decrypt the message sent from Bob. Meanwhile, Bob needs to prove “Bob is Bob” because if an attacker Eve utilized Alice's public key to encrypt forged messages and sent them to Alice, then Alice could get wrong messages from Eve (in this case, Eve is the fake “Bob”). In this consideration, Bob needs to utilize his own private key to form a digital signature based on the one side hash function of the original message and send it to Alice. Because the digital signature includes the information of Bob's private key, then the message of path 2 can be proven to have come from the real Bob. To verify whether the message from path 1 does not contain counterfeit information, the message decrypted from Bob's side is also passed through Hash Function to compare with the decrypted messages from path 2. Using this approach, the information security and the authentication could be guaranteed.

Here, an example of how to generate the asymmetric random keys in RSA cryptosystem is shown. Supposing two simple prime number 5 and 11 are chosen for the RSA-based cryptosystem, a public key and private key of (55, 3) and (55, 27) could be generated. If Bob wants to deliver a message of “26” then a cyphertext of 16 is delivered to Alice. Alice could use the private key to decode. If an attacker Eve wants to decode the decoding process, she should find the D number and calculate the 26²⁷ mod 55, which is 160059109085386090080713531498405298176. Thus, with only small prime numbers 5 and 11 already being able to generate this large product, choosing much larger prime numbers could further improve the robustness of the cryptosystem. For more information on the RSA algorithm, key generation process and calculation details can be found in Supplementary Materials Note I and Table S1. Based on the above analysis and calculation, the robustness of the cryptosystem is based on both the difficulty of solving discrete logarithms and the randomness of the memristive entropy source. Moreover, the intrinsic stochasticity of the spatiotemporal variability of the memristor provides temporally-changing, physically-unpredictable and machine-learning-attack-resilient primitives for random key generation, which also prevents the forward secrecy (FS) issue in the cryptography field. The FS is to ensure that the session key could not be compromised if the historical secrets have been compromised. Even if the attacker gets the session key one time, the previous and forward key generation is still highly secure. Moreover, the currently proposed memristive TRNG could be employed against the side channel attacks. As shown in Fig. S7, the physical power consumption cannot derive any useful information of the random numbers.

In this work, memristive spatiotemporal variation based on a graphene CBRAM has been comprehensively studied. Towards cryptographic key generation, three levels of analysis have been carried out: physical level, AI level, and random number level. On the physical-level analysis, a KMC simulation has been carried out to demonstrate the revolution of the filament during forming, Set and Reset, and the DDV and CCV are visualised. On the AI level, PCA, LSTM and SHAP value techniques have been utilized to investigate whether the intrinsic memristive stochastic parameters are truly random. The random-number-level analysis is based on Hamming distances calculation, which verifies the probability, stability, reproducibility and uniqueness of those random number based on the memristive primitives. Apart from the Set voltage distributions of the DDV and CCV, the RTN-based TTV has been investigated by kernel density map and time constant, which indicates equal chances of the upper and lower current levels. Merging the spatiotemporal variation including DDV, CCV and TTV as hardware primitive enhances the robustness of the cryptographic key generation. For practical application, Rivest-Shamir-Adleman (RSA) cryptosystem based on the memristive spatiotemporal variation has been demonstrated. As proven in the current study, the memristive spatiotemporal variability provides true randomness and also improves the robustness of the cryptosystem.

Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.chip.2023.100040.

Acknowledgments This study was supported by grants from National Natural Science Foundation of China (62174008), Beijing Municipal Education Commission (KZ202110005001), the Ministry of Science and Technology, Taiwan, China (MOST 111-2119-M-492- 002-MBK, MOST 111-2622-8-182-001-TS1, MOST 109-2221-E-182-013-MY3, and MOST 110- 2221-E-182-043-MY3), and the Chang Gung Memorial Hospital (CORPD2J0073).

Declaration of Competing Interest The authors declare no competing interests.