Mathematical Biosciences 315 (2019) 108218
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Termite population size estimation based on termite tunnel patterns using a convolutional neural network
T
Jeong-Kweon Seoa, Seongbok Baikb, Sang-Hee Leec,d,
⁎
a
Department of Mathematics, College of Natural Sciences, Chonnam National University, 77, Yongbong-ro, Buk-gu, Gwangju 61186, Republic of Korea KT Network Laboratory, 1689-70 Yuseong St., Yuseong, Daejeon 34047, Republic of Korea c Division of Integrated Mathematics, National Institute for Mathematical Sciences, Daejeon, 34047, Republic of Korea d Center for Convergent Research of Emerging Virus Infection, Korea Research Institute of Chemical Technology, Daejeon 34114, Republic of Korea b
ARTICLE INFO
ABSTRACT
Keywords: Termite tunnel pattern Convolutional neural network Population size estimation Agent-based model
Subterranean termites live in large colonies under the ground where they build an elaborate network of tunnels for foraging. In this study, we explored how the termite population size can be estimated using partial information on tunnel patterns. To achieve this, we used an agent-based model to create tunnel patterns that were characterized by three variables: the number of simulated termites (N), passing probability of two termites encountering one another (P), and distance that termites move soil particles (D). To explore whether the N value could be estimated using a partial termite tunnel pattern, we generated four tunnel pattern groups by partially obscuring different areas in an image of a complete tunnel pattern, where: (1) the outer area of the tunnel pattern was obscured (I-pattern); (2) half of the tunnel pattern was obscured (H-pattern); (3) the inner region of the tunnel pattern was obscured (O-pattern); and (4) I- and O-patterns (IO-pattern) were combined. For each group, 80% of the tunnel patterns were used to train a convolutional neural network while the remaining 20% were used for estimating the N value. Estimation results showed that the N estimates for IO-patterns were the most accurate, followed by I-, H-, and O-patterns. This indicates that termite population size can be estimated based on tunnel information near the center of a colony. We briefly discuss the advantages and disadvantages of this method for estimating termite population size.
1. Introduction Subterranean termites (specifically Coptotermes or some Rhinotermitidae) construct underground tunnel networks for foraging [1–3]. They contribute to ecosystem stability by recycling rotten wood and other cellulose-based materials. In contrast, they invade and destroy wooden buildings constructed at ground level, which results in economic damage. The scale of the damage depends on how many termites live below the ground where the damage occurred, with more individuals probably causing greater damage than fewer individuals. To date, several methods have been used to estimate population size of termite colonies formed under the ground. In the case of mound- and tree-nesting termite species, a direct counting method has been used [4]. Labeling insects with radioactive isotopes has been a popular insect-marking method. However, these methods can not be used directly to estimate population size of subterranean termites because they mainly perform the foraging activities under concrete slabs and inside
buildings [5]. In the case of subterranean termite species, the mark-release-recapture (MRR) method [6,7] and multiple capture-recapture (MCR) method [8] have been widely used. The MCR method repeats the process of capturing, marking, releasing, and recapturing termites to investigate how many of the marked individuals were caught. In many cases, the MCR method is more accurate than the MRR method, as the latter method only performs the process once. The merits of MRR and MCR are that they can be applied directly to subterranean termite species and do not adversely affect their colonies [9–11]. On the other hand, these two methods have disadvantages in that a large error may occur in termite population estimates depending on the state of the surrounding environment. [12]. These methods also include mathematical assumptions that contribute to accuracy errors. One of the assumptions is that the capture probabilities of marked termites are identical across populations at separate points in time [13–15]. In real termite colony systems, ecological factors such as differential trap
⁎ Corresponding author at: National Institute for Mathematical Sciences, KT Daeduk, 2 Research Center, 463-1, Jeonmin-dong, Yusung-gu, Daejeon, Republic of Korea. E-mail address:
[email protected] (S.-H. Lee).
https://doi.org/10.1016/j.mbs.2019.108218 Received 24 December 2018; Received in revised form 17 June 2019; Accepted 18 June 2019 Available online 19 June 2019 0025-5564/ © 2019 Published by Elsevier Inc.
Mathematical Biosciences 315 (2019) 108218
J.-K. Seo, et al.
exposure, and the heterogeneity of the soil substrate, are likely to violate some assumptions [16]. Additionally, many marked termites often die because of various environmental factors, which results in an error in the estimation of the population size [17]. To overcome the disadvantages of MRR and MCR, Seo et al. [18] published a method for estimating the number of termites present in a tunnel network based on information on tunnel patterns. In their study, the authors simulated termite tunnel patterns using an agent-based model (ABM). Next, they studied the tunnel patterns using a convolutional neural network (CNN) and attempted to estimate the N values for the new tunnel patterns. However, this method is very difficult to apply to actual termite colonies, necessitating the need for more information on the structure of an entire termite tunnel pattern. Considering the huge size of termite colonies, directly digging out the soil was not an option. Thus, in this study, we investigated whether the N value could be better estimated by only using a part, instead of an entire tunnel pattern, to improve on the idea proposed by Seo et al. [18]. In addition, we briefly discuss the advantages of the proposed methods and issues to be overcome in future studies.
narrow tunnel, they show one of three types of response behavior. In the first response, the two proceed in the advancing direction without physical interaction. This behavior was reported in a study by Sim and Lee [23]. In the second type of response, one of the two termites turns and walks in the opposite direction. In the third type of behavior, one of the two termites begins to excavate a new tunnel to the right or left side wall. In this ABM, the probability (P) that each behavior type would be selected was set at 0.9 for passing, 0.05 for turning, and 0.05 for tunneling [23]. In preliminary studies, we checked that the termite tunneling patterns, generated when the P value was within the range of (0.7–0.9), did not change significantly. 2.5. Analysis Initially, we generated termite tunnel patterns with 100 × 100 resolution using the ABM. However, these patterns were resized to 30 × 30 resolution for application in the CNN using different model parameter values. We trained tunnel patterns with 30 × 30 resolution for application in a CNN. For resizing, we used the resize function provided by MATLAB software (Mathworks Inc., Natick, MA). The CNN consisted of two parts, feature extraction and classification, in which the weighted parameters were calculated during the entire training procedure [24]. The feature extraction part consisted of a sub-sampling (pooling) layer and a convolution layer (see Fig. 1). Each layer was represented by a two-dimensional map of weights linked together. The convolutional training process that used the shared weights of the input data was as follows:
2. Materials and methods 2.1. Model description The ABM used in this study was also used by Seo et al. [18] and Song and Lee [19]. A 100 × 100 grid space was used to simulate termite movement. One simulated termite occupies two grid spaces. One grid represents the head and the other represents the body. We set all simulated termites to overlap at the center of the simulation space at the start of the simulation (t = 0). And we created a circular empty area with a radius of 10 grids in the center of the space. This area was initially a place for the simulated termites to tunnel into and collect soil moved from the end of a tunnel during tunneling activity.
H 1 H 1
uij =
Xi + p, j + q × hpq, q=1 p=1
where Xij denotes the pixel position of (i, j) in the tunnel pattern image to be input for learning, and the value represents the image intensity of the corresponding pixel. Here, the tunnel has a value of 1 and the empty space has a value of 0. h denotes a convolution filter of size H × H (=3). In the process of training the convolutional layers, the stride value was set to two. After the convolutional computation, a pooling process was implemented. Among several pooling methods, we used the following max pooling method:
2.2. Traveling behavior of simulated termites Simulated termites could either walk, carry sand partices from the end of a tunnel, and create new tunnels. Actual temrites change their direction regularly to the left and right as they move along the tunnel. This can be understood as a behavior to obtain geometry information of the tunnel wall [21]. To mimic these behaviors, the simulated termites adopted a rule that allowed them to move randomly to the left and right when moving straight ahead. In addition, we included a rule in the ABM to eliminate the loop-back structure of tunnels that are not present in real termite tunnels. This rule prevented the simulated termites from creating tunnels with loops. Experiments in homogeneous sand reveal that termites do not usually make loops. However, a loop can be formed in the process of tunnel-tunnel intersection [22].
uij =
max Xpq
(p, q) Pij
where Pij indicates a set of pixel indices at (i, j) locations in the predefined pooling receptive field. To ensure the pooling layer or convolution layer to generate the appropriate output size, we used zeropadding for the pooling output. The length of the zero-padding unit was given by ((filter size - 1)) / 2. For the convolution layer and the pooling layer, a rectified linear unit (ReLu) function was used for the non-
2.3. Transport behavior of simulated termites Once the simulated termites pick up the sand particles from the tip of the tunnel, they turn against the direction of the body and moved toward the center of the nest. The termites move 20 steps and then drop the sand particle on the site. Then they walk back toward the end of the tunnel. In this model, the steps were used as the variable D and the value was fixed as the value of 20. In preliminary studies, when the D value was between 15 and 35, it was confirmed that there was almost no effect on the tunnel pattern. Through this process, termite tunnels grow over time [20]. With this ABM, we allowed the soil moved by the termites to be deposited in the same grid site, to prevent the tunnel from being clogged by soil. 2.4. Interactions between encountering termites Fig. 1. Structure of the convolutional neural network (CNN) used in this study.
When two actual termites move forward facing each other in a 2
Mathematical Biosciences 315 (2019) 108218
J.-K. Seo, et al.
Table 1 Structure of CNN as a list of layers and activations. The filter size and the number of filters are represented as n × n and m, respectively (the detailed structure is shown in Fig. 1). Layer
(n × n) × m
Activation
Conv. Conv. Pooling Conv. Conv. Pooling Conv. Fully-Conn. Dropout Fully-Conn.
(3 × 3) × 10 (3 × 3) × 100 Kernel size: (2 × 2) (3 × 3) × 210 (3 × 3) × 441 Kernel size: (2 × 2) (3 × 3) × 882 Output units: 256 Dropout rate: 0.5 Output units: 7
ReLu ReLu Max strides: 2 ReLu ReLu Max strides: 2 ReLu ReLu Softmax
Fig. 3. Typical tunnel patterns generated by simulated termites. N is the number of simulated termites.
saturating activation function f(x). ReLu not only has no effect on the receptive field of the convolutional layer but also has a positive role in improving the nonlinear nature of the decision function and the overall network. The final layer of the CNN structure was connected to a multilayer perceptron neural network with an appropriate dropout ratio (=0.5). This dropout function prevented the overfitting by cutting off the specific connections to the neural network. Following this forward passing calculation process in each training epoch, the calculation of the backpropagation process was performed in the reverse direction of the CNN architecture. This updated the minimal gradient changes of the weighted parameter values in the kernels and classification neural network. Details of the architecture are shown in Table 1. We set the threshold for the softmax function output to 0.5. If the Nth node of the softmax output exceeded the threshold value, then the N value of the tunnel pattern image was calculated. If the values of all nodes were less than the threshold value, then the image estimate was ignored. Implementing the numerical experiment that applied CNN generated 7000 tunnel patterns with grid sizes of 30 × 30 for three different parameter values: the number of simulated termites (N), the passing probability of two encountering termites (P), and the distance moved by termites to deposit soil parcels during tunneling activity (D). We resized tunnel patterns with 100 × 100 resolution to reduce tunnel pattern learning time. Fig. 2 shows the original tunnel pattern image, revealing that this ABM can produce a pattern similar to the actual
termite tunnel pattern [18]. Four types of tunnel pattern images were created by covering parts of the tunnel patterns, where: (1) the outer area of the tunnel pattern was obscured and only the tunnel structure of the inner area was shown (I-pattern), (2) half of the tunnel pattern was obscured (H-pattern), (3) the inner region of the tunnel pattern was obscured (O-pattern), and (4) I- and O-patterns (IO-pattern) were combined (see Fig. 2). The four tunnel pattern groups were created to investigate whether N values could be estimated from partly obscured tunnel patterns. We trained CNN using 80% of the patterns (800 for N = 3, 4, …, 9, respectively) and estimated N values for the remaining 20% patterns (200 for N = 3, 4, …, 9, respectively). 3. Results Fig. 3 shows the simulated termite tunnel patterns for N = 3, 4… 9. As N increased, the tunnel patterns became more complicated. This tendency was evident in the range of N = 3 to 6, but not for N = 7, 8, and 9. In particular, the complexity near the center of the tunnel increased significantly as N increased. This was largely due to two reasons: a larger number of simulated termites increases the number of primary tunnels formed at the beginning of tunneling, and greater termite numbers in the tunnels increase the likelihood of them encountering one other, which also increases the frequency in which two termites block each other's way. Blocking can result in new tunnels, which contribute to the complexity of the entire tunnel pattern. The estimation accuracy (ε) is defined as the number of tunnels tested and the number of tunnels in which the N value was estimated correctly. The error rate generated when measuring actual termite population size led us to assume that the correct values of N, N + 1, and N−1 had been used for the simulation. In this study, we generated four tunnel pattern groups by obscuring parts of the tunnel patterns in the complete tunnel pattern image (I-, H-, O- and IO-pattern; Fig. 2). Fig. 4 shows the ratio of learning to validation accuracy when CNN learnt each tunnel pattern group. Learning was optimized when the ratio value was minimized. For I- and H- patterns, the learning rate of CNN quickly reached the optimized state, where it was slightly slower than for IO-pattern. Additionally, CNN took the longest time for O-pattern. After the initial optimization state, there were several similar states. This indicates that CNN learning was not as good for O-patterns when compared with other patterns. Table 2 shows the results of estimating the N values for O-pattern, H-pattern, I-, and IO-pattern. For O-pattern, CNN demonstrated that when N = 3 and 4, the ε was approximately 72% and 73%, respectively. However, for N = 4 and 5, the values of ε were lower than 50%. When N > 5, the values of ε were more than 50%. These results suggest
Fig. 2. Entire and partially obscured termite tunnel patterns. I-pattern corresponds to the central region of the entire pattern, H-pattern is half of the entire pattern, O-pattern is the obscured central area of the entire pattern, and IOpattern corresponds to the central region and partial outer region of the entire pattern. 3
Mathematical Biosciences 315 (2019) 108218
J.-K. Seo, et al.
Fig. 4. Training and validation accuracy of the convolutional neural network (CNN). Arrows indicate maximum accuracies.
for simulating termite tunnel patterns only considers walkers, not soldiers. However, since Reticulitermes and Coptotermes are known to have 2% and 10% of soldiers in their colonies, respectively, we can estimate total population size in these termite colonies. The estimation method proposed in this study is based on two models. One is a convolutional neural network (CNN), a type of deep learning process, while the other is an agent-based model (ABM). The CNN was used to extract tunnel pattern characteristics and provided an advantage in that we did not need to define exact tunnel pattern features. The ABM generated tunnel pattern image data for CNN learning. However, these two models had some limitations. The CNN limitations are related to the classification being dependent on the number of hidden layers, the selection of the pooling method, and resolution of the input image. Therefore, it was necessary to select an optimal CNN structure for maximizing classification performance. In this study, we did not investigate the optimal condition because it showed high accuracy for estimating the N value even under rough conditions. As for the ABM limitations, the simulated tunnel patterns obtained from this model may be different from actual tunnel patterns created by termites in the field. This is due to some ABM parameter values being based on experimental data obtained from a homogeneous sand substrate [21,25]. The soil wherein termites actually live is typically heterogeneous due to the diversity of soil properties. Termites have very different foraging behaviors depending on soil moisture concentration and soil particle density. Furthermore, these behaviors differ depending on the termite species. If soil conditions are not suitable for termite activity, then tunneling behavior is reduced and tunnel lengths are shortened [26]. Termite behavior is also strongly influenced by the roughness of the tunnel's base surface [27]. In addition to heterogeneity effects, a lack of experimental data was a weakness in the algorithm for determining simulated termite behavior. Simulated termites are supposed to drop soil particles from the tunnel tip at a distance of D, which was a constant value. However, actual termites seem to drop soil after they have moved a variety of distances, rather than a specific distance. Our ongoing experiments have found that these divergent distances seem to follow a specific distribution function, which will be examined in future studies. Moreover, we fixed the P and D values at 0.9 and 20, respectively. In
Table 2 Estimation accuracy of the trained convolutional neural network (CNN) for O-, H-, I-, and IO-patterns with N = 3, 4… 9. N used in simulations
Estimation accuracy (%) O-patterns H-patterns
I-patterns
IO-patterns
3 4 5 6 7 8 9
71.8 73.2 47.0 41.8 50.2 66.1 59.1
85.7 84.6 80.1 66.8 63.9 77.4 70.3
86.0 85.9 76.2 76.9 64.6 72.6 68.8
76.7 80.0 60.6 61.2 59.3 68.7 55.0
that O-patterns contain very little characteristic tunnel information that can be distinguished according to N values. For H-pattern, the value of ε had a relatively high value for all cases except N = 9 when compared with the learning of O-pattern. Interestingly, for both O- and H-patterns, the value of ε was very high for N = 3 and 4. This indicates that simpler tunnel patterns provide more characteristic information than complex patterns, which varies depending on the N value. For I-patterns, the ε values were higher for all N values than that for H-patterns. 4. Discussion and conclusion The MRR and MCR methods have been mainly used to estimate termite population size. However, disadvantages of these methods include the amount of time needed to obtain an estimate and error values greatly fluctuating due to different environmental conditions. In our previous study [18], we proposed a new approach to overcome these shortcomings. This method was based on overall tunnel pattern information. However, the large sizes of the tunnels made it difficult to obtain the entire tunnel pattern information. The present study is a follow-up on Seo et al. [18], which explored whether it was possible to estimate population size using partial tunnel pattern information. To be precise, this study does not estimate termite population size in a colony but estimates the number of tunnel excavators. This is because the ABM 4
Mathematical Biosciences 315 (2019) 108218
J.-K. Seo, et al.
preliminary studies, we checked that the termite tunneling patterns, generated when the P and D values were within the range of (0.7–0.9) and (15–35), respectively, did not change significantly. However, depending on the environmental conditions, these values may be outside this range; thus, further research is required to determine how tunnel pattern varies with changes in these variable values. The fundamental problem with the ABM model is the limited simulation space size and number of excavators. In other words, it is questionable whether the ABM can describe the actual termite tunnel patterns more precisely when the number of excavators and space size are increased. Thus said, we do not believe that the ABM can generate a more elaborate tunnel pattern even if this condition is met because the variable values used in the model in this study were based on experiments with approximately 10,000 termites. In this study, N = 6 corresponded to 10,000 actual termites. Our earlier work [18] addressed this problem. A larger termite population within the tunnel will result in more complex individual-individual interactions, which is likely to cause the emergence of greater tunnel patterns. We recommend that experiments with various termite population sizes is essential to solve the problems associated with space and population size. If the new variables obtained from these experiments are taken into account in the current model, then the model used in this study may be a more useful tool for estimating termite population size. Although the models used in this study have several drawbacks, our study is considered valuable not only because it presents a new approach to estimating termite population size, but it also provides the ABM that can be used to study termite tunnel patterns. In particular, this model could be used to improve the baiting methods used for termite pest control. Thus far, termite researchers have tried to remove termite colonies by equally spacing bait stations within a certain distance around the point where termite damage has occurred. However, this method has several disadvantages such as environmental pollution and excessive costs. Using the termite tunnel pattern information obtained from the model used in this study, we can explore optimal baitstation installation methods that will permit bait to spread more quickly and effectively to all termites in the nest.
References [1] E.C. King, W.T. Spink, Foraging galleries of the Formosan subterranean termite, Coptotermes formosanus, in Louisiana, Ann. Entomol. Am. 62 (1969) 536. [2] G. Henderson, H. Fei, Comparison of native subterranean termite and Formosan subterranean termite: biology, ecology, and methods of control, Forest Products Society Conference, Kissimmee, Florida, 2002, pp. 11–13. [3] N.-Y. Su, M. Tamashiro, J.R. Yates, M.I. Haverty, Foraging behavior of the Formosan subterranean termite (Isoptera: Rhinotermitidae), Environ. Entomol. 13 (1984) 1466. [4] G. Josens, K. Soki, Relation between termite numbers and the size of their mounds, Insectes Soc. 57 (2010) 303. [5] N.-Y. Su, R.H. Scheffrahn, Foraging population and territory of the formosan subterranean termite (Isoptera: rhinotermitidae) in an urban environment, Sociobiology 14 (1988) 353. [6] K.H. Pollock, J.D. Nichols, C. Brownie, J.E. Hines, Statistical inference for capture–recapture experiments, Wildlife Monographs 107 (1990) 97. [7] B.T. Forschler, M.L. Townsend, Mark-release-recapture estimates of reticulitermes spp. (Isoptera: rhinotermitidae) colony foraging populations from Georgia, USA Environ. Entomol. 25 (1996) 952. [8] W. Hayes, R. Carter, Population monitoring, in: A. Alberts (Ed.), West Indian Iguanas: Status Survey and Conservation Action Plan, IUCN, Gland, Switzerland, 2000, pp. 79–85. [9] T. Nobre, L. Nunes, D.E. Bignell, Estimation of foraging territories of Reticulitermes grassei through mark-release-recapture, Entomol. Exp. Appl. 123 (2007) 119. [10] A. Arab, A.M. Costa-Leonardo, F.E. Casarin, A.C. Guaraldo, R.C. Chaves, Foraging activity and demographic patterns of two termite species (Isoptera: rhinotermitidae) living in urban landscapes in southeastern Brazil, Eur. J. Entomol. 102 (2005) 691. [11] K.S. Brown, G.H. Broussard, B.M. Kard, A.L. Smith, M.P. Smith, Colony characterization of Reticulitermes flavipes (Isoptera: Rhinotermitidae) on a native tallgrass prairie, Am. Midl. Nat. 159 (2008) 21. [12] M.W.J Crosland, N.-Y. Su, Mark-recapture without estimating population sizes: a tool to evaluate termite baits, Bull. Entomol. Res. 96 (2006) 99. [13] H. Puche, N.-Y. Su, Estimating population density of the Formosan subterranean termite, Coptotermes formosanus(Isoptera: Rhinotermitidae) using the effective sampling area of in-ground monitoring stations, Bull. Entomol. Res. 94 (2004) 47. [14] S.-H. Lee, N.-Y. Su, Estimating the size of a closed population from a single markrecapture regime by using a diffusion probability of marked individuals, Appl. Math. Comput. 182 (2006) 403. [15] N.-Y. Su, S.-H. Lee, Estimating the population size and colony boundary of subterranean termites by using the density functions of directionally-averaged capture probability, J. Econ. Entomol. 101 (2008) 592. [16] R.M. Dorazio, J.A. Royle, Mixture models for estimating the size of a closed population when capture rates vary among individuals, Biometrics 59 (2003) 351. [17] B.L. Thorne, E. Russek-Cohen, B.T. Forschler, N.L. Breisch, J.F.A. Traniello, Evaluation of mark-release-recapture methods for estimating forager population size of subterranean termite (Isoptera: Rhinotermitiae) colonies, Environ. Entomol. 25 (1996) 938. [18] J.-K. Seo, S. Baik, S.-H. Lee, Simulating the architecture of a termite incipient nest using a convolutional neural network, Ecol. Inform. 44 (2018) 94. [19] H.-S. Song, S.-H. Lee, Simulation study to determine why only some termites are active during tunneling activity, Entomol. Sci. 21 (2018) 185. [20] H.-F. Li, N.-Y. Su, Sand displacement during tunnel excavation by the formosan subterranean termite (Isoptera: Rhinotermitidae), Ann. Entomol. Soc. Am. 10 (2008) 456. [21] S.W. Sim, S.-H. Kang, S.-H. Lee, Using hidden Markov models to characterize termite traveling behavior in tunnels with different curvatures, Behav. Processes 111 (2015) 101. [22] P. Bardunias, N.-Y. Su, Dead reckoning in tunnel propagation of the formosan subterranean termite (Isoptera: rhinotermitidae), Ann. Entomol. Soc. Am. 102 (2009) 158. [23] S.W. Sim, S.-H. Lee, Measurement of the time required for termites to pass each other in tunnels of different curvatures, J. Insect Sci. 20 (2013) 550. [24] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (2015) 436. [25] N.-Y. Su, B.M. Stith, H. Puche, P. Bardunias, Characterization of tunneling geometry of subterranean termites (Isoptera: rhinotermitidae) by computer simulation, Sociobiology 44 (2004) 471. [26] N.-Y. Su, H. Puche, Tunneling activity of subterranean termites (Isoptera: rhinotermitidae) in sand with moisture gradients, J. Econ. Entomol. 96 (2003) 88. [27] S.-H. Lee, Traveling behavior of a termite in tunnels with different curvatures and based surface roughness, J. Asia Pac. Entomol. 21 (2018) 258.
Declaration of Competing Interest None. Acknowledgments This work was supported by the National Institute for Mathematical Sciences(NIMS) and also suppored by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (2017R1E1A1A03070059). Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mbs.2019.108218.
5