Comparison of two commercial DNA extraction kits for the analysis of nasopharyngeal bacterial communities

Marcos Pérez-Losada; Keith A. Crandall; Robert J. Freishtat; Marcos Pérez-Losada; Keith A. Crandall; Robert J. Freishtat

doi:10.3934/microbiol.2016.2.108

AIMS Microbiology

2016, Volume 2, Issue 2: 108-119. doi: 10.3934/microbiol.2016.2.108

Previous Article Next Article

Communication Topical Sections

Comparison of two commercial DNA extraction kits for the analysis of nasopharyngeal bacterial communities

1.
Computational Biology Institute, George Washington University, Ashburn, VA 20147, USA
2.
Division of Emergency Medicine, Children’s National Medical Center, Washington, DC 20010, USA
3.
CIBIO-InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, Universidade do Porto, Campus Agrário de Vairão, Vairão 4485-661, Portugal

Received: 31 March 2016 Accepted: 28 April 2016 Published: 29 April 2016

Characterization of microbial communities via next-generation sequencing (NGS) requires an extraction ofmicrobial DNA. Methodological differences in DNA extraction protocols may bias results and complicate inter-study comparisons. Here we compare the effect of two commonly used commercial kits (Norgen and Qiagen)for the extraction of total DNA on estimatingnasopharyngeal microbiome diversity. The nasopharynxis a reservoir for pathogens associated with respiratory illnesses and a key player in understandingairway microbial dynamics.
Total DNA from nasal washes corresponding to 30 asthmatic children was extracted using theQiagenQIAamp DNA and NorgenRNA/DNA Purification kits and analyzed via IlluminaMiSeq16S rRNA V4 ampliconsequencing. The Norgen samples included more sequence reads and OTUs per sample than the Qiagen samples, but OTU counts per sample varied proportionallybetween groups (r = 0.732).Microbial profiles varied slightly between sample pairs, but alpha- and beta-diversity indices (PCoAand clustering) showed highsimilarity between Norgen and Qiagenmicrobiomes. Moreover, no significant differences in community structure (PERMANOVA and adonis tests) and taxa proportions (Kruskal-Wallis test) were observed betweenkits. Finally, aProcrustes analysis also showed low dissimilarity (M² = 0.173; P< 0.001) between the PCoAs of the two DNA extraction kits.
Contrary to what has been observed in previous studies comparing DNA extraction methods, our 16S NGS analysis of nasopharyngeal washes did not reveal significant differences in community composition or structure between kits. Our findingssuggest congruence between column-based chromatography kits and supportthe comparison of microbiomeprofilesacross nasopharyngeal metataxonomic studies.

Keywords:

Citation: Marcos Pérez-Losada, Keith A. Crandall, Robert J. Freishtat. Comparison of two commercial DNA extraction kits for the analysis of nasopharyngeal bacterial communities[J]. AIMS Microbiology, 2016, 2(2): 108-119. doi: 10.3934/microbiol.2016.2.108

Related Papers:

[1]	Nurul Akhma Zakaria, A.A. Kutty, M.A. Mahazar, Marina Zainal Abidin . Arsenic acute toxicity assessment on select freshwater organism species in Malaysia. AIMS Environmental Science, 2016, 3(4): 804-814. doi: 10.3934/environsci.2016.4.804
[2]	Patricia Morcillo, Maria Angeles Esteban, Alberto Cuesta . Mercury and its toxic effects on fish. AIMS Environmental Science, 2017, 4(3): 386-402. doi: 10.3934/environsci.2017.3.386
[3]	Antonio Zanin, Ivonez Xavier de Almeida, Francieli Pacassa, Fabricia Silva da Rosa, Paulo Afonso . Maturity level of environmental management in the pulp and paper supply chain. AIMS Environmental Science, 2021, 8(6): 580-596. doi: 10.3934/environsci.2021037
[4]	Djordje Vilimanovic, Gangadhar Andaluri, Robert Hannah, Rominder Suri, A. Ronald MacGillivray . Occurrence and aquatic toxicity of contaminants of emerging concern (CECs) in tributaries of an urbanized section of the Delaware River Watershed. AIMS Environmental Science, 2020, 7(4): 302-319. doi: 10.3934/environsci.2020019
[5]	Zeynab Rezazadeh Salteh, Saeed Fazayeli, Saeid Jafarzadeh Ghoushchi . Evaluation and prioritization of barriers to the implementation of the eco-regenerative supply chains using fuzzy ZE-numbers framework in group decision-making. AIMS Environmental Science, 2024, 11(4): 516-550. doi: 10.3934/environsci.2024026
[6]	Francesco Teodori . Health physics calculation framework for environmental impact assessment of radiological contamination. AIMS Environmental Science, 2021, 8(4): 403-420. doi: 10.3934/environsci.2021026
[7]	P. U. Singare, J.P. Shirodkar . Persistent and carcinogenic polycyclic aromatic hydrocarbons in the north-western coastal marine environment of India. AIMS Environmental Science, 2021, 8(2): 169-189. doi: 10.3934/environsci.2021012
[8]	Melanie Voigt, Indra Bartels, Anna Nickisch-Hartfiel, Martin Jaeger . Elimination of macrolides in water bodies using photochemical oxidation. AIMS Environmental Science, 2018, 5(5): 372-388. doi: 10.3934/environsci.2018.5.372
[9]	Fernanda P. Maciel, Joshua M. Peschel . A GIS-based tool for bioaccumulation risk analysis and its application to study polychlorinated biphenyls in the Great Lakes. AIMS Environmental Science, 2018, 5(1): 1-23. doi: 10.3934/environsci.2018.1.1
[10]	Jamshid Ali . Environmental resilience: transition to regenerative supply chain management. AIMS Environmental Science, 2024, 11(2): 107-128. doi: 10.3934/environsci.2024007

Abstract

1. Introduction

With the continuous development and progress of computing technology over the past decade, deep learning has taken a big step on the road of evolution. Deep learning can be found in many fields, such as imaging and natural language. Scholars have begun to apply deep learning to solve some complex partial differential equations (PDEs), which include PDEs with high order derivatives ^[1], high-dimensional PDEs ^[2,3], subdiffusion problems with noisy data ^[4] and so on. Based on the deep learning method, Raissi et al. ^[5] proposed a novel algorithm called the physics informed neural network (PINN), which has made excellent achievements for solving forward and inverse PDEs. It integrates physical information described by PDEs into a neural network. In recent years, the PINNs algorithm has attracted extensive attention. To solve forward and inverse problems of integro-differential equations (IDEs), Yuan et al. ^[6] proposed auxiliary physics informed neural network (A-PINN). Lin and Chen ^[7] designed a two-stage physics informed neural network for approximating localized wave solutions, which introduces the measurement of conserved quantities in stage two. Yang et al. ^[8] developed Bayesian physics informed neural networks (B-PINNs) that takes Bayesian neural network and Hamiltonian Monte Carlo or the variational inference as a priori and posteriori estimators, respectively. Scholars have also presented other variant algorithms of PINN, such as RPINNs ^[9] and $ hp $-VPINNs ^[10]. PINN has also achieved good performance for solving physical problems, including high-speed flows ^[11] and heat transfer problems ^[12]. In addition to integer-order differential equations, authors have studied the application of PINN in solving fractional differential equations such as fractional advection-diffusion equations (see Pang et al. ^[13] for fractional physics informed neural network (fPINNs)), high dimensional fractional PDEs (see Guo et al. ^[14] for Monte Carlo physics-informed neural networks (MC-PINNs)) and fractional water wave models (see Liu et al. ^[15] for time difference PINN).

In recent decades, fractional differential equations (FDEs) have been concerned and studied in many fields, such as image denoising ^[16] and physics ^[17,18,19]. The reason why fractional differential equations have attracted wide attention is that they can more clearly describe complex physical phenomena. As an indispensable part of fractional problems, distributed-order differential equations are difficult to solve due to the complexity of distributed-order operators. To solve the time multi-term and distributed-order fractional sub-diffusion equations, Gao et al. ^[20] proposed a second order numerical difference formula. Jian et al. ^[21] derived a fast second-order implicit difference scheme of time distributed-order and Riesz space fractional diffusion-wave equations and analyzed the unconditional stability and second-order convergence. Li et al. ^[22] applied the mid-point quadrature rule with finite volume method to approximate the distributed-order equation. For the nonlinear distributed-order sub-diffusion model ^[23], the distributed-order derivative and the spatial direction were approximated by the FBN-$ \theta $ formula with a second-order composite numerical integral formula and the $ H^1 $-Galerkin mixed finite element method, respectively. In ^[24], Guo et al. adopted the Legendre-Galerkin spectral method for solving 2D distributed-order space-time reaction-diffusion equations. For the two-dimensional Riesz space distributed-order equation, Zhang et al. ^[25] used Gauss quadrature to calculate the distributed-order derivative and applied an alternating direction implicit (ADI) Galerkin-Legendre spectral scheme to approximate the spatial direction. For the distributed-order fourth-order sub-diffusion equation, Ran and Zhang ^[26] developed new compact difference schemes and proved their stability and convergence. In ^[27,28], authors developed spectral methods for the distributed-order time fractional fourth-order PDEs.

As we know, distributed-order fractional PDEs can be regarded as the limiting case of multi-term fractional PDEs ^[29]. Moreover, Diethelm and Ford ^[30] have observed that small changes in the order of a fractional PDE lead to only slight changes in the final solution, which gives initial support to the employed numerical integration method. In view of this, we develop the FBN-$ \theta $ ^[31,32] with a second-order composite numerical integral formula combined with a multi-output neural network for solving 1D and 2D nonlinear time distributed-order models. Based on the idea of using a single output neural network combined with the discrete scheme of fractional models ^[13], we also use a single output neural network combined with a time discrete scheme to solve the nonlinear time distributed-order models. However, the accuracy of the prediction solution calculated by the single output neural network scheme is low and the training progress takes a lot of time. Therefore, we introduce a multi-output neural network to obtain the numerical solution of the time discrete scheme. Compared with the single output neural network scheme, the proposed multi-output neural network scheme has two main advantages as follows:

● Saving more computing time. The multi-output neural network scheme makes the sampling domain of the collocation points from the spatiotemporal domain to the spatial domain, which decreases the number of the training dataset and thus reduces the training time.

● Improving the accuracy of predicted solution. Due to the discrete scheme of the distributed-order derivative, the $ n $-th output item of the multi-output neural network will be constrained by the previous $ n-1 $ output items.

The remainder of this article is as follows: In Section 2, we show what the components of neural network are and how to construct a neural network. In Section 3, we give the lemmas used to approximate the distributed-order derivative and the process of building the loss function. In Section 4, we provide some numerical results to confirm the capability of our proposed method. Finally, we make some conclusions in Section 5.

2. Neural network

In face of different objectives in various fields, scholars have developed many different types of neural networks, such as feed-forward neural network (FNN) ^[6], recurrent neural network (RNN) ^[33] and convolutional neural network (CNN) ^[34]. The FNN considered in this article can effectively solve most PDEs. Input layer, hidden layer and output layer are three indispensable components of FNN, which can be given, respectively, by

$\begin{equation} \notag \begin{aligned} \text{input layer:}\; \; \; \Phi_0(x) & = x\in\mathbb{R}^{d_{in}}, \\ \text{hidden layers:}\; \; \; \Phi_k(x) & = \sigma(\boldsymbol W_k\Phi_{k-1}(x)+\boldsymbol b_k)\in\mathbb{R}^{\lambda_k}, \; 1\leq k\leq K-1, \\ \text{output layer:}\; \; \Phi_K(x) & = \boldsymbol W_K\Phi_{K-1}(x)+\boldsymbol b_K\in\mathbb{R}^{d_{out}}. \end{aligned} \end{equation}$ $

$ \boldsymbol W_k\in\mathbb{R}^{\lambda_k\times\lambda_{k-1}} $ and $ \boldsymbol b_k\in\mathbb{R}^{\lambda_k} $ represent the weight matrix and the bias vector in the $ k $th layer, respectively. We define $ \delta = \{\boldsymbol W_k, \boldsymbol b_k\}_{1\leq k\leq K} $, which is the trainable parameters of FNN. $ \lambda_k $ represents the number of neurons included in the $ k $th layer. $ \sigma $ is a nonlinear activation function. In this article, the hyperbolic tangent function ^[3,6] is selected as the activation function. There are many other functions that can be considered as activation functions, such as the rectified linear unit (ReLU) $ \sigma(x) = \max\{x, 0\} $ ^[4] and the logistic sigmoid $ \sigma(x) = \frac{1}{1+e^{-x}} $ ^[35].

3. Methodology

3.1. Problem setup

In this article, we consider a nonlinear distributed-order model with the following general form:

$\begin{equation} \begin{split} \left\{ \begin{aligned} & D^w_tu+\mathcal{N}(u) = f(\boldsymbol x, t), \; (\boldsymbol x, t) \in \Omega \times J, \\ & u(\boldsymbol x, 0) = u_0(\boldsymbol x), \boldsymbol x\in\bar{\Omega}, \end{aligned} \right. \end{split} \end{equation}$ $

(3.1)

where $ \Omega \subset \mathbb{R}^d(d\leq2) $ and $ J = (0, T] $. $ \mathcal{N}[\cdot] $ is a nonlinear differential operator. $ D^w_tu $ represents the distributed-order derivative and has the following definition:

$\begin{equation} \begin{split} D^w_tu(\boldsymbol x, t) = \int^1_0 \omega(\alpha){}^C_0D^{\alpha}_tu(\boldsymbol x, t)d\alpha, \end{split} \end{equation}$ $

(3.2)

where $ \omega(\alpha) \geq 0 $, $ \int^1_0\omega(\alpha)d\alpha = c_0 > 0 $ and $ {}^C_0D^{\alpha}_tu(\boldsymbol x, t) $ is the Caputo fractional derivative expressed by

$\begin{equation} {}^C_0D^{\alpha}_tu(\boldsymbol x, t) = \left\{ \begin{aligned} & \frac{1}{\Gamma(1-\alpha)}\int^t_0\frac{ u_{\eta}(\boldsymbol x, \eta)}{(t-\eta)^{-\alpha}}d\eta, \; 0 < \alpha < 1, \\ & u_t(\boldsymbol x, t), \alpha = 1. \end{aligned} \right. \end{equation}$ $

(3.3)

The specific boundary condition is determined by the practical problem.

3.2. Some lemmas

For simplicity, choosing a mesh size $ \Delta \alpha = \frac{1}{2I} $, we denote the nodes on the interval $ [0, 1] $ with coordinates $ \alpha_i = i\Delta \alpha $ for $ i = 0, 1, 2, \cdots, 2I $. The time interval $ [0, T] $ is divided as the uniform mesh with the grid points $ t_n = n\Delta t (n = 0, 1, 2, \cdots, N) $ and $ \Delta t = T/N $ is the time step size. We denote $ v^n \approx u^n = u(\boldsymbol x, t_n) $, $ u^{n+\frac{1}{2}}_t = \frac{u^{n+1} - u^n}{\Delta t} + O(\Delta t^2) $ and $ u^{n+\frac{1}{2}} : = \frac{u^{n+1}+u^n}{2} $. $ v^n $ is defined as the approximation solution of the time discrete scheme. The following lemmas are introduced to construct the numerical discrete formula of (3.1):

Lemma 3.1. (See ^[23]) Supposing $ \omega(\alpha) \in C^2[0, 1] $, we can get

$\begin{equation} \int^1_0\omega(\alpha)d\alpha = \Delta \alpha\mathop \sum \limits_{k = 0}^{2I} c_i\omega(\alpha_i)-\frac{\Delta \alpha^2}{12}\omega^{(2)}(\gamma), \; \gamma \in (0, 1), \end{equation}$ $

(3.4)

where

$\begin{equation} c_i = \left\{ \begin{aligned} & \frac{1}{2}, \; i = 0, 2I, \\ & 1, \; \mathit{\text{otherwise}}. \end{aligned} \right. \end{equation}$ $

(3.5)

Lemma 3.2. From ^[23,31,32], the discrete formula of the Caputo fractional derivative (3.3) can be obtained by

$\begin{equation} {}^C_0D^{\alpha}_tu(\boldsymbol x, t_{n+\frac{1}{2}}) = \frac{{}^C_0D^{\alpha}_tu^{n+1} + {}^C_0D^{\alpha}_tu^n}{2} + O(\Delta t^2) = \Delta t^{-\alpha}\mathop \sum \limits_{s = 0}^{n + 1} \widetilde{\kappa}^{(\alpha)}_{n+1-s}u^s + O(\Delta t^2), \end{equation}$ $

(3.6)

where

$\begin{equation} \widetilde{\kappa}^{(\alpha)}_{n+1-s} = \left\{ \begin{aligned} & \frac{\kappa^{(\alpha)}_{0}}{2}, \; s = n+1, \\ & \frac{\kappa^{(\alpha)}_{n-s} + \kappa^{(\alpha)}_{n+1-s}}{2}, \; otherwise. \end{aligned} \right. \end{equation}$ $

(3.7)

The parameters $ \kappa^{(\alpha)}_{i}(i = 0, 1, \cdots, n+1) $ that are the coefficients of FBN-$ \theta $ $ (\theta \in [-\frac{1}{2}, 1]) $ can be given by

$\begin{equation} \kappa^{(\alpha)}_{i} = \left\{ \begin{aligned} & 2^{-\alpha}(1+\alpha\theta)(3-2\theta)^{\alpha}, \; i = 0, \\ & \frac{\phi_0\kappa^{(\alpha)}_{0}}{\psi_0}, \; i = 1, \\ & \frac{1}{2\psi_0}[(\phi_0-\psi_1)\kappa^{(\alpha)}_{1} + \phi_1\kappa^{(\alpha)}_{0}], \; i = 2, \\ & \frac{1}{i\psi_0}\sum^{3}_{j = 1}[\phi_{j-1}-(i-j)\psi_j]\kappa^{(\alpha)}_{i-j}, \; i\geq3, \end{aligned} \right. \end{equation}$ $

(3.8)

where

$\begin{equation} \phi_i = \left\{ \begin{aligned} & 2\alpha(\theta-1)(\alpha\theta+1)+\alpha\theta(\theta-\frac{3}{2}), \; i = 0, \\ & -\alpha(2\theta^2-3\alpha\theta+4\alpha\theta^2-1), \; i = 1, \\ & -\alpha\theta(\frac{1}{2}-\theta+\alpha-2\alpha\theta), \; i = 2, \\ \end{aligned} \right. \end{equation}$ $

(3.9)

and

$\begin{equation} \psi_i = \left\{ \begin{aligned} & \frac{1}{2}(3-2\theta)(1+\alpha\theta), \; i = 0, \\ & -\frac{\alpha\theta}{2}(3-2\theta)-2(1-\theta)(\alpha\theta+1), \; i = 1, \\ & -\frac{1}{2}(2\theta-1)(\alpha\theta+1)-2\alpha\theta(\theta-1), \; i = 2, \\ & -\frac{1}{2}\alpha\theta(1-2\theta). \end{aligned} \right. \end{equation}$ $

(3.10)

Lemma 3.3. (See ^[23]) The distributed-order term $ D^{\omega}_tu $ at $ t = t_{n+\frac{1}{2}} $ can be calculated by the following formula:

$\begin{equation} D^{\omega}_tu(\boldsymbol x, t_{n+\frac{1}{2}}) = \frac{D^{\omega}_tu^{n+1} + D^{\omega}_tu^{n}}{2}+O(\Delta t^2) = \frac{1}{2}\mathop \sum \limits_{s = 0}^{n + 1} \beta^{n+1}_su^s+O(\Delta t^2+\Delta\alpha^2), \end{equation}$ $

(3.11)

where

$\begin{equation} \beta^{n+1}_s = \left\{ \begin{aligned} & \widehat{\kappa}_{n-s}+\widehat{\kappa}_{n+1-s}, \; 0\leq s < n+1, \\ & \widehat{\kappa}_0, \; s = n+1, \\ \end{aligned} \right. \end{equation}$ $

(3.12)

and

$\begin{equation} \widehat{\kappa}_{n-s} = \mathop \sum \limits_{i = 0}^{2I} \frac{\varphi_i}{\Delta t^{\alpha_i}}\kappa^{(\alpha_i)}_{n-s}, \; \varphi_i = \Delta\alpha\omega(\alpha_i)c_i. \end{equation}$ $

(3.13)

3.3. The loss function

Based on the above lemmas, the discrete scheme of the distributed-order model (3.1) at $ t = t_{n+\frac{1}{2}}(n = 0, 1, 2, \cdots, N-1) $ can be expressed by the following equality:

$\begin{equation} \frac{1}{2}\mathop \sum \limits_{s = 0}^{n + 1} \beta^{n+1}_sv^s + \mathcal{N}(v^{n+\frac{1}{2}}) = f^{n+\frac{1}{2}}. \end{equation}$ $

(3.14)

Then we can obtain the system of equations as follows:

$\begin{equation} \left\{ \begin{aligned} & \frac{1}{2}\sum^{1}_{s = 0}\beta^{1}_sv^s + \mathcal{N}(v^{\frac{1}{2}}) = f^{\frac{1}{2}}, \\ & \frac{1}{2}\sum^{2}_{s = 0}\beta^{2}_sv^s + \mathcal{N}(v^{1+\frac{1}{2}}) = f^{1+\frac{1}{2}}, \\ & \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \vdots \\ & \frac{1}{2}\sum^{N}_{s = 0}\beta^{N}_sv^s + \mathcal{N}(v^{N-\frac{1}{2}}) = f^{N-\frac{1}{2}}. \end{aligned} \right. \end{equation}$ $

(3.15)

The system of Eq (3.15) can be rewritten as the following matrix form:

$\begin{equation} \boldsymbol v(\boldsymbol x)M + \boldsymbol{\mathcal{N}}(\boldsymbol v(\boldsymbol x)) + v^0(\boldsymbol x)\boldsymbol\rho_0 = \boldsymbol f(\boldsymbol x), \end{equation}$ $

(3.16)

where the symbols $ \boldsymbol\rho_0 $, $ \boldsymbol v(\boldsymbol x) $, $ \boldsymbol f(\boldsymbol x) $ and $ \boldsymbol{\mathcal{N}}(\boldsymbol v(\boldsymbol x)) $ are vectors, which are given, respectively, by

$\begin{equation} \notag \begin{aligned} & \boldsymbol\rho_0 = [\frac{1}{2}\beta^1_0, \frac{1}{2}\beta^2_0, \frac{1}{2}\beta^3_0, \cdots, \frac{1}{2}\beta^N_0], \\ \\ & \boldsymbol v(\boldsymbol x) = [v^1(\boldsymbol x), v^2(\boldsymbol x), \cdots, v^N(\boldsymbol x)], \\ \\ & \boldsymbol f(\boldsymbol x) = [\boldsymbol f^{\frac{1}{2}}(\boldsymbol x), \boldsymbol f^{1+\frac{1}{2}}(\boldsymbol x), \cdots, \boldsymbol f^{N-\frac{1}{2}}(\boldsymbol x)], \\ \\ & \boldsymbol{\mathcal{N}}(\boldsymbol v(\boldsymbol x)) = [\mathcal{N}(v^{\frac{1}{2}}(\boldsymbol x)), \mathcal{N}(v^{1+\frac{1}{2}}(\boldsymbol x)), \cdots, \mathcal{N}(v^{N-\frac{1}{2}}(\boldsymbol x))]. \end{aligned} \end{equation}$ $

The symbol $ M $ is an $ N\times N $ matrix that has the following definition:

$\begin{equation} \notag \begin{aligned} M = \begin{pmatrix} \frac{1}{2}\beta^1_1 & \frac{1}{2}\beta^2_1 &\cdots &\frac{1}{2}\beta^N_1 \\ 0 & \frac{1}{2}\beta^2_2 & \cdots &\frac{1}{2}\beta^N_2 \\ \vdots& \vdots& &\vdots\\ 0 & 0& \cdots & \frac{1}{2}\beta^N_N \end{pmatrix}. \end{aligned} \end{equation}$ $

Now, we introduce a multi-output neural network $ \boldsymbol v(\boldsymbol x; \delta) = [v^1(\boldsymbol x; \delta), v^2(\boldsymbol x; \delta), \cdots, v^N(\boldsymbol x; \delta)] $ into Eq (3.16), which takes $ \boldsymbol x $ as an input and is used to approximate time discrete solutions $ \boldsymbol v(\boldsymbol x) = [v^1(\boldsymbol x), v^2(\boldsymbol x), \cdots, v^N(\boldsymbol x)] $. This will result in a multi-output PINN $ \boldsymbol \ell(\boldsymbol x) = [\ell^1(\boldsymbol x), \ell^2(\boldsymbol x), \cdots, \ell^N(\boldsymbol x)] $:

$\begin{equation} \boldsymbol \ell(\boldsymbol x) = \boldsymbol v(\boldsymbol x; \delta)M + \boldsymbol{\mathcal{N}}(\boldsymbol v(\boldsymbol x; \delta)) + v^0(\boldsymbol x)\boldsymbol\rho_0 - \boldsymbol f(\boldsymbol x), \end{equation}$ $

(3.17)

where $ \ell^{n+1}(\boldsymbol x) $ is denoted as residual error of the discrete scheme (3.14), which is given by

$\begin{equation} \ell^{n+1}(\boldsymbol x) = \frac{1}{2}\mathop \sum \limits_{s = 1}^{n + 1} \beta^{n+1}_sv^s(\boldsymbol x; \delta) + \mathcal{N}(v^{n+\frac{1}{2}}(\boldsymbol x; \delta)) + \frac{1}{2}\beta^{n+1}_0v^0(\boldsymbol x)-f^{n+\frac{1}{2}}(\boldsymbol x), \; n = 0, 1, 2, \cdots, N-1. \end{equation}$ $

(3.18)

The loss function is constructed in the form of mean square error. Combined with the boundary condition loss, the total loss function can be expressed by the following formula:

$\begin{equation} MSE_{total} = MSE_{\ell}+MSE_b, \end{equation}$ $

(3.19)

where

$\begin{equation} MSE_{\ell} = \frac{1}{N\times N_x}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_x}} |\ell^j(\boldsymbol x^i_{\ell})|^2, \end{equation}$ $

(3.20)

and boundary condition loss

$\begin{equation} MSE_b = \frac{1}{N\times N_b}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_b}} |v^j(\boldsymbol x^i_b; \delta)-u^j(\boldsymbol x^i_b)|^2. \end{equation}$ $

(3.21)

Here, $ \{\boldsymbol x^i_{\ell}\}^{N_x}_{i = 1} $ corresponds to the collocation points on the space domain $ \Omega $ and $ \{\boldsymbol x^i_{b}\}^{N_b}_{i = 1} $ denotes the boundary training data. The schematic diagram of using the multi-output neural network scheme to solve nonlinear time distributed-order models is shown in Figure 1.

Figure 1. A multi-output neural network framework to solve nonlinear time distributed-order models, where $ MSE_* $ represents the loss function.

DownLoad: Full-Size Img PowerPoint

4. Algorithm implementation

In this section, we consider two nonlinear time distributed-order equations to verify the feasibility and effectiveness of our proposed method. The performance is evaluated by calculating the relative $ L^2 $ error between the predicted and exact solutions. The definition of relative $ L^2 $ error is given by

$\begin{equation} ||u-v||_{L^2} = \frac{\sqrt{\mathop \sum \nolimits_{j = 1}^N \mathop \sum \nolimits_{i = 1}^{{N_x}} |u^j(\boldsymbol x^i)-v^j(\boldsymbol x^i)|^2}}{\sqrt{\sum^N_{j = 1}\mathop \sum \nolimits_{i = 1}^{{N_x}} |u^j(\boldsymbol x^i)|^2}}. \end{equation}$ $

(4.1)

Table 1 indicates which optimizer is selected for each example to minimize the loss function.

Table 1. The hyperparameters configured in numerical examples.

Example	Optimizer	Learning rate	Iterations
$ 1 $	Adam + L-BFGS	0.001	20000
$ 2 $	Adam + L-BFGS	0.001	20000
$ 3 $	Adam + L-BFGS	0.001	20000
$ 4 $	L-BFGS	-	-

| Show Table

DownLoad: CSV

We use Python to code our algorithms and all codes run on a Lenovo laptop with AMD R7-6800H CPU @ 3.20 GHz and 16.0GB RAM.

4.1. The distributed-order sub-diffusion model

Here, we solve the following distributed-order sub-diffusion model:

$\begin{equation} u_t+D^{\omega}_tu-\Delta u-\Delta u_t+G(u) = f(\boldsymbol x, t), \; (\boldsymbol x, t)\in\Omega\times J, \end{equation}$ $

(4.2)

with boundary condition

$\begin{equation} u(\boldsymbol x, t) = 0, \; \boldsymbol x\in\partial\Omega, \; t\in\bar{J}, \end{equation}$ $

(4.3)

and initial condition

$\begin{equation} u(\boldsymbol x, 0) = u_0(\boldsymbol x), \; \boldsymbol x\in\bar{\Omega}, \end{equation}$ $

(4.4)

where the nonlinear term $ G(u) = u^2 $. The symbol $ \Delta $ is the Laplace operator. Based on Eqs (3.14)–(3.21), the loss function $ MSE_{total} $ can be obtained by

$\begin{equation} \notag MSE_{total} = MSE_{\ell}+MSE_b, \end{equation}$ $

where

$\begin{equation} \notag \begin{aligned} & MSE_{\ell} = \frac{1}{N\times N_x}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_x}} |\ell^j(\boldsymbol x^i_{\ell})|^2, \\ & MSE_b = \frac{1}{N\times N_b}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_b}} |v^j(\boldsymbol x^i_b; \delta)-0|^2. \end{aligned} \end{equation}$ $

Example 1.

For this example, we set space domain $ \Omega = (0, 1) $ and time interval $ J = (0, \frac{1}{2}] $. The training set consists of the boundary points and $ N_x = 200 $ collocation points randomly selected in the space domain $ \Omega $. Choosing $ \omega(\alpha) = \Gamma(3-\alpha) $ and the source term

$\begin{equation} \notag f(x, t) = 2t\sin(2\pi x)+\frac{\Gamma(3) t(t-1)}{\ln t}\sin(2\pi x)+4t^2\pi^2\sin(2\pi x)+8t\pi^2\sin(2\pi x)+(t^2\sin(2\pi x))^2, \end{equation}$ $

then we can obtain the exact solution $ u(x, t) = t^2\sin(2\pi x) $.

To evaluate the performance of our proposed method, the exact solution and the predicted solution solved by a multi-output neural network that consists of $ 6 $ hidden layers with 40 neurons in each hidden layer are showed in Figure 2.

Figure 2. Example 1: the exact solution and predicted solution with $ \Delta \alpha = \frac{1}{500} $, $ \theta = 1 $ and $ N = 20 $ at $ t = 0.5 $.

DownLoad: Full-Size Img PowerPoint

The influence of different network structures on our proposed method to solve Example 1 is presented in Table 2. The accuracy of the predicted solution fluctuates with different network architectures, and shows a fluctuating growth with the increase of the number of hidden layers. Based on the three network architectures, Figure 3 shows the behavior performance of the proposed method with gradual decrease of the time step size. With the increase of the number of grid points in the time interval, we observe that the behavior of relative $ L^2 $ error generally presents an upward trend for fixed network architecture and expanding the depth of the hidden layer can effectively improve the accuracy of the predicted solution.

Table 2. The relative $ L^2 $ error between the predicted solution with parameters $ N = 20 $, $ \Delta\alpha = \frac{1}{500} $, $ \theta = 1 $ and exact solution for different numbers of hidden layers and different numbers of neurons per layer.

	20	30	40	50	60
2	3.583482e-02	7.850782e-02	5.238403e-02	1.143600e-01	3.761048e-02
4	3.815002e-02	5.144618e-02	3.622479e-02	2.317038e-02	3.721057e-02
6	2.433509e-02	1.216386e-01	2.638224e-02	2.837553e-02	2.682958e-02

| Show Table

DownLoad: CSV

Figure 3. Example 1: the variation trend of relative $ L^2 $ error between the predicted solution with $ \theta = 1 $, $ \Delta\alpha = \frac{1}{500} $ and the exact solution.

DownLoad: Full-Size Img PowerPoint

Numerical results calculated by the single output neural network and multi-output neural network schemes are presented in Table 3, where we select $ 200 $ collocation points in the given spatial domain by random sampling method and set $ N = 10 $, $ \theta = 1 $ and $ \Delta\alpha = \frac{1}{500} $. It is easy to find that the accuracy of the predicted solution calculated by the multi-output neural network scheme is higher than that solved by the single output neural network scheme and replacing single output neural network with multi-output neural network can save a lot of computing time.

Table 3. The relative $ L^2 $ error and computing time given by the single output and multi-output neural network schemes.

Neural network	Layers	Neurons	Relative $ L^2 $ error	CPU time (s)
multi-output	4	20	1.948281e-02	32.80
	4	40	1.858120e-02	44.45
	6	20	9.724722e-03	51.40
	6	40	1.024869e-02	65.88
single output	4	20	8.741336e-02	661.99
	4	40	2.702814e-02	742.17
	6	20	2.937448e-02	748.31
	6	40	2.278774e-02	840.09

| Show Table

DownLoad: CSV

Example 2.

In this numerical example, considering the space domain $ \Omega = (0, 1)\times(0, 1) $, the time interval $ J = (0, \frac{1}{2}] $, $ \omega(\alpha) = \Gamma(3-\alpha) $ and the exact solution $ u(x, y, t) = t^2\sin(2\pi x)\sin(2\pi y) $, the source term can be given by

$\begin{equation} \notag \begin{aligned} f(x, y, t) = & 2t\sin(2\pi x)\sin(2\pi y) + \frac{\Gamma(3)t(t-1)}{\ln t}\sin(2\pi x)\sin(2\pi y) \\ & + 8t^2\pi^2\sin(2\pi x)\sin(2\pi y) + 16t\pi^2\sin(2\pi x)\sin(2\pi y) + (t^2\sin(2\pi x)\sin(2\pi y))^2. \end{aligned} \end{equation}$ $

The training dataset is shown in Figure 5 and the collocation and boundary points are selected by random sampling method.

To better illustrate the behavior of the predicted solution, Figure 4 portrays the contour plot of the exact and predicted solutions, where the training set consists of $ 961 $ collocation points and $ 124 $ boundary points selected by the equidistant uniform sampling method and the network architecture consists of 12 hidden layers with 60 neurons per layer.

Figure 4. Example 2: the contour plot of the exact (a) and predicted (b) solutions with $ N = 10 $, $ \theta = 1 $ and $ \Delta\alpha = \frac{1}{500} $ at $ t = 0.5 $.

DownLoad: Full-Size Img PowerPoint

Figure 5. Example 2: distribution of the collocation points (a) sampled in the domain $ \Omega $ and the boundary training dataset (b).

DownLoad: Full-Size Img PowerPoint

Table 4 shows the impact of depth and width of the network on the accuracy of the predicted solution. In Figure 6, we present the behavior of how relative $ L^2 $ error changes with respect to different grid points $ N $. Combined with Table 4 and Figure 6, we observe that increasing the number of hidden layers or neurons has a positive effect on reducing relative $ L^2 $ error in general.

Table 4. The relative $ L^2 $ error between the predicted solution with parameters $ N = 20 $, $ \theta = 1 $, $ \Delta\alpha = \frac{1}{500} $ and the exact solution for different number of hidden layers and neurons per layer.

	20	30	40	50	60
4	1.155984e-01	1.101123e-01	7.228454e-02	7.787707e-02	4.492646e-02
6	6.273340e-02	6.088659e-02	9.802474e-02	6.536637e-02	6.245142e-02
8	6.163051e-02	7.424834e-02	6.844676e-02	5.500003e-02	5.169303e-02

| Show Table

DownLoad: CSV

Figure 6. Example 2: the variation trend of relative $ L^2 $ error between the predicted solution with $ \theta = 1 $, $ \Delta\alpha = \frac{1}{500} $ and the exact solution.

DownLoad: Full-Size Img PowerPoint

The results shown in Table 5 reveal the performance of the single output neural network and multi-output neural network schemes, where we select $ 40 $ boundary points and $ 200 $ collocation points in the given spatial domain by random sampling method and set $ N = 10 $, $ \theta = 1 $ and $ \Delta\alpha = \frac{1}{500} $. One can see that using multi-output neural network effectively improves the precision and reduces the computing time.

Table 5. The relative $ L^2 $ error and computing time given by the single output and multi-output neural network schemes.

Neural network	Layers	Neurons	Relative $ L^2 $ error	CPU time (s)
multi-output	4	20	1.080369e-01	52.84
	4	40	7.472573e-02	66.05
	6	20	6.654115e-02	72.71
	6	40	7.937601e-02	93.77
single output	4	20	3.443264e-01	731.22
	4	40	7.054495e-01	728.11
	6	20	2.636401e-01	772.08
	6	40	1.410780e-01	892.94

| Show Table

DownLoad: CSV

4.2. The distributed-order fourth-order sub-diffusion model

Further, we consider the following distributed-order fourth-order sub-diffusion model:

$\begin{equation} u_t+D^{\omega}_tu-\Delta u+\Delta^2u+G(u) = f(\boldsymbol x, t), \; (\boldsymbol x, t)\in\Omega\times J, \end{equation}$ $

(4.5)

with boundary condition

$\begin{equation} \notag u(\boldsymbol x, t) = \Delta u(\boldsymbol x, t) = 0, \; \boldsymbol x\in\partial\Omega, \; t\in\bar{J}, \end{equation}$ $

and initial condition

$\begin{equation} \notag u(\boldsymbol x, 0) = u_0(\boldsymbol x), \boldsymbol x\in\bar{\Omega}, \end{equation}$ $

where the nonlinear term $ G(x) = u^2 $.

Similarly, the corresponding loss function $ MSE_{total} $ can be calculated by

$\begin{equation} \notag MSE_{total} = MSE_{\ell}+MSE_b, \end{equation}$ $

where

$\begin{equation} \notag \begin{aligned} MSE_{\ell} = & \frac{1}{N\times N_x}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_x}} |\ell^j(\boldsymbol x^i_{\ell})|^2, \\ MSE_b = &\frac{1}{N\times N_b}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_b}} |v^j(\boldsymbol x^i_b; \delta)-0|^2 \\ & + \frac{1}{N\times N_b}\mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{i = 1}^{{N_b}} |\Delta v^j(\boldsymbol x^i_b; \delta)-0|^2. \end{aligned} \end{equation}$ $

Example 3.

Here, we define the space-time domain $ \Omega\times J = (0, 1)\times(0, \frac{1}{2}] $. Considering $ \omega(\alpha) = \Gamma(3-\alpha) $ and the source term

$\begin{equation} f(x, t) = 2t\sin(\pi x)+\frac{\Gamma(3)t(t-1)}{\ln t}\sin(\pi x)+t^2\pi^2\sin(\pi x)+t^2\pi^4\sin(\pi x)+(t^2\sin(\pi x))^2, \end{equation}$ $

(4.6)

the exact solution can be given by $ u(x, t) = t^2\sin(\pi x) $. Similar to Example $ 1 $, we also randomly sample $ 200 $ collocation points in the space domain $ \Omega $.

In order to conveniently observe the capability of our proposed method, Figure 7 shows the change in the trajectory of the predicted and exact solutions with respect to the space point $ x $, where the parameters are set as $ N = 20 $, $ \Delta\alpha = \frac{1}{500} $, $ \theta = 1 $ and the network consists of $ 6 $ hidden layers with $ 50 $ neurons per layer. Figure 8 portrays the trajectory of relative $ L^2 $ error with three different network architectures. Table 6 shows the impact of expanding depth or width of the network on the accuracy of the predicted solutions. Based on Figure 8 and Table 6, it is easy to observe that increasing the number of hidden layer plays a positive role in improving the accuracy of the predicted solutions.

Figure 7. Example 3: the predicted and exact solutions at $ t = 0.5 $.

DownLoad: Full-Size Img PowerPoint

Figure 8. Example 3: the relative $ L^2 $ error between the predicted solution with $ \theta = 1 $, $ \Delta \alpha = \frac{1}{500} $ and exact solution for different number of grid points.

DownLoad: Full-Size Img PowerPoint

Table 6. The relative $ L^2 $ error between the predicted solution with parameters $ N = 20 $, $ \theta = 1 $, $ \Delta\alpha = \frac{1}{500} $ and the exact solution for different number of hidden layers and neurons per layer.

	20	30	40	50	60
2	8.994465e-03	8.835675e-03	7.976196e-03	1.121301e-02	1.284193e-02
4	6.698531e-03	6.623791e-03	1.012124e-02	7.129064e-03	6.760145e-03
6	5.238653e-03	3.178300e-03	4.705057e-03	7.221552e-03	5.002096e-03

| Show Table

DownLoad: CSV

The relative $ L^2 $ error and CPU time obtained by the multi-output neural network and single output neural network schemes are presented in Table 7, where we select $ 200 $ collocation points in the given spatial domain by random sampling method and set $ N = 10 $, $ \theta = 1 $ and $ \Delta\alpha = \frac{1}{500} $. The error of the proposed multi-output neural network scheme is smaller than that of the single output neural network scheme. For this 1D system, the multi-output neural network scheme is more efficient than the single output neural network scheme.

Table 7. The relative $ L^2 $ error and computing time given by the single output and multi-output neural network schemes.

Neural network	Layers	Neurons	Relative $ L^2 $ error	CPU time (s)
multi-output	4	20	8.027377e-03	334.07
	4	40	8.656712e-03	366.52
	6	20	8.315435e-03	557.78
	6	40	6.788274e-03	601.19
single output	4	20	2.412753e-02	994.62
	4	40	2.365781e-02	1317.16
	6	20	2.428662e-02	1344.76
	6	40	2.139562e-02	1795.75

| Show Table

DownLoad: CSV

Example 4.

Now we take space domain $ \Omega = (0, 1)\times(0, 1) $ and time interval $ J = (0, \frac{1}{2}] $. Let $ \omega(\alpha) = \Gamma(3-\alpha) $ and the exact solution $ u(x, y, t) = t^2\sin(\pi x)\sin(\pi y) $. Then we arrive at the source term

$\begin{equation} \begin{aligned} f(x, y, t) = & 2t\sin(\pi x)\sin(\pi y) + \frac{\Gamma(3)t(t-1)}{\ln t}\sin(\pi x)\sin(\pi y)+2t^2\pi^2\sin(\pi x)\sin(\pi y)\\ & +4t^2\pi^4\sin(\pi x)\sin(\pi y)+(t^2\sin(\pi x)\sin(\pi y))^2. \end{aligned} \end{equation}$ $

(4.7)

Here, we apply the training data set shown in Figure 5.

In order to more intuitively demonstrate the feasibility of our proposed method for solving this 2D system, Figure 9 shows the contour plot of $ u $ and $ |u-v| $, where the training set consists of $ 900 $ collocation points and $ 120 $ boundary points selected by the equidistant uniform sampling method and the network is composed of 6 hidden layers with 20 neurons per layer.

Figure 9. Example 4: the contour plot of $ u $ and $ |u-v| $ with $ \theta = 1 $, $ \Delta \alpha = \frac{1}{500} $ and $ N = 40 $ at $ t = 0.5 $.

DownLoad: Full-Size Img PowerPoint

From the behavior of relative $ L^2 $ error in Figure 10, one can see that the accuracy of the predicted solutions with the fixed network first presents a increasing trend and then gradually presents a downward trend. This is because the approximation ability of neural network reaches saturation with the increase of grid points $ N $. Table 8 shows the relative $ L^2 $ error calculated by different network architectures. On the whole, the relative $ L^2 $ error slightly decreases with expanding depth of the network, while it first slightly decreases and then increases with expanding width of the network. To show the precision and efficiency of the multi-output neural network scheme for this 2D system, the relative $ L^2 $ error and CPU time obtained by the multi-output neural network and single output neural network schemes are shown in Table 9, where we select $ 40 $ boundary points and $ 200 $ collocation points in the given spatial domain by random sampling method and set $ N = 10 $, $ \theta = 1 $ and $ \Delta\alpha = \frac{1}{500} $. It illustrates that the multi-output neural network scheme is more accurate and efficient than the single output neural network scheme.

Figure 10. Example 4: the relative $ L^2 $ error between the predicted solution with $ \theta = 1 $, $ \Delta \alpha = \frac{1}{500} $ and the exact solution for different number of grid points.

DownLoad: Full-Size Img PowerPoint

Table 8. The relative $ L^2 $ error between the predicted solution with parameters $ N = 40 $, $ \theta = 1 $, $ \Delta\alpha = \frac{1}{500} $ and the exact solution for different number of hidden layers and neurons per layer.

}	20	30	40	50	60
2	7.593424e-02	4.556132e-02	2.926755e-02	5.850593e-02	4.956921e-02
4	3.716823e-02	2.498269e-02	3.388639e-02	4.073820e-02	5.282780e-02
6	2.690082e-02	3.039336e-02	2.769315e-02	4.331023e-02	6.456580e-02

| Show Table

DownLoad: CSV

Table 9. The relative $ L^2 $ error and computing time given by the single output and multi-output neural network schemes.

Neural network	Layers	Neurons	Relative $ L^2 $ error	CPU time (s)
multi-output	4	20	4.549939e-02	345.03
	4	40	5.847813e-02	339.46
	6	20	5.680785e-02	466.02
	6	40	6.110534e-02	528.10
single output	4	20	1.455785e-01	902.99
	4	40	1.330922e-01	1309.05
	6	20	1.496581e-01	1469.85
	6	40	1.322566e-01	1990.96

| Show Table

DownLoad: CSV

5. Conclusions

In this article, a multi-output physics informed neural network combined with the Crank-Nicolson scheme including the FBN-$ \theta $ method and the composite numerical integral formula was constructed to solve 1D and 2D nonlinear time distributed-order models. The calculation process is described in detail. Numerical experiments are provided to prove the effectiveness and feasibility of our algorithm. Compared with the results calculated by a single output neural network combined with the FBN-$ \theta $ method and the Crank-Nicolson scheme, one can clearly see that the proposed multi-output neural network scheme is more efficient and accurate. Moreover, some numerical methods, such as finite difference or finite element method, need to linearize the nonlinear term, which will give rise to extra costs. The process of linearization can be directly omitted by PINN. Further work will investigate the application of the proposed methodology in high-dimensional problems and practical problems ^{[36,37,38,39,40]}.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors would like to thank the editor and all the anonymous referees for their valuable comments, which greatly improved the presentation of the article. This work is supported by the National Natural Science Foundation of China (12061053, 12161063), Natural Science Foundation of Inner Mongolia (2021MS01018), Young innovative talents project of Grassland Talents Project, Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region (NMGIRT2413, NMGIRT2207), and 2023 Postgraduate Research Innovation Project of Inner Mongolia (S20231026Z).

Conflict of interest

The authors declare that they have no conflict of interest.

References

[1]	Ding T, Schloss PD (2014) Dynamics and associations of microbial community types across the human body. Nature 509: 357–360. doi: 10.1038/nature13178
[2]	Human Microbiome Project C (2012) Structure, function and diversity of the healthy human microbiome. Nature 486: 207–214. doi: 10.1038/nature11234
[3]	Human Microbiome Project C (2012) A framework for human microbiome research. Nature 486: 215–221. doi: 10.1038/nature11209
[4]	Integrative HMPRNC (2014) The Integrative Human Microbiome Project: dynamic analysis of microbiome-host omics profiles during periods of human health and disease. Cell Host Microbe 16: 276–289. doi: 10.1016/j.chom.2014.08.014
[5]	Marchesi JR, Ravel J (2015) The vocabulary of microbiome research: a proposal. Microbiome 3: 31. doi: 10.1186/s40168-015-0094-5
[6]	Kuczynski J, Lauber CL, Walters WA, et al. (2012) Experimental and analytical tools for studying the human microbiome. Nat Rev Genet 13: 47–58.
[7]	Althani AA, Marei HE, Hamdi WS, et al. (2015) Human Microbiome and its Association With Health and Diseases. J Cell Physiol 231: 1688–1694.
[8]	Martin R, Miquel S, Langella P, et al. (2014) The role of metagenomics in understanding the human microbiome in health and disease. Virulence 5: 413–423. doi: 10.4161/viru.27864
[9]	Cox MJ, Cookson WO, Moffatt MF (2013) Sequencing the human microbiome in health and disease. Hum Mol Genet 22: R88–94. doi: 10.1093/hmg/ddt398
[10]	Brooks JP, Edwards DJ, Harwich MD, Jr., et al. (2015) The truth about metagenomics: quantifying and counteracting bias in 16S rRNA studies. BMC Microbiol 15: 66. doi: 10.1186/s12866-015-0351-6
[11]	Abusleme L, Hong BY, Dupuy AK, et al. (2014) Influence of DNA extraction on oral microbial profiles obtained via 16S rRNA gene sequencing. J Oral Microbiol 6.
[12]	Wu GD, Lewis JD, Hoffmann C, et al. (2010) Sampling and pyrosequencing methods for characterizing bacterial communities in the human gut using 16S sequence tags. BMC Microbiol 10: 206. doi: 10.1186/1471-2180-10-206
[13]	Momozawa Y, Deffontaine V, Louis E, et al. (2011) Characterization of bacteria in biopsies of colon and stools by high throughput sequencing of the V2 region of bacterial 16S rRNA gene in human. PLOS ONE 6: e16952. doi: 10.1371/journal.pone.0016952
[14]	Lazarevic V, Gaia N, Girard M, et al. (2013) Comparison of DNA extraction methods in analysis of salivary bacterial communities. PLOS ONE 8: e67699. doi: 10.1371/journal.pone.0067699
[15]	Willner D, Daly J, Whiley D, et al. (2012) Comparison of DNA extraction methods for microbial community profiling with an application to pediatric bronchoalveolar lavage samples. PLOS ONE 7: e34605. doi: 10.1371/journal.pone.0034605
[16]	Bogaert D, De Groot R, Hermans PW (2004) Streptococcus pneumoniae colonisation: the key to pneumococcal disease. Lancet Infect Dis 4: 144–154. doi: 10.1016/S1473-3099(04)00938-7
[17]	Garcia-Rodriguez JA, Fresnadillo Martinez MJ (2002) Dynamics of nasopharyngeal colonization by potential respiratory pathogens. J Antimicrob Chemother 50 Suppl S2: 59–73.
[18]	Biesbroek G, Tsivtsivadze E, Sanders EA, et al. (2014) Early respiratory microbiota composition determines bacterial succession patterns and respiratory health in children. Am J Respir Crit Care Med 190: 1283–1292. doi: 10.1164/rccm.201407-1240OC
[19]	Teo SM, Mok D, Pham K, et al. (2015) The infant nasopharyngeal microbiome impacts severity of lower respiratory infection and risk of asthma development. Cell Host Microbe 17: 704–715. doi: 10.1016/j.chom.2015.03.008
[20]	Feazel LM, Santorico SA, Robertson CE, et al. (2015) Effects of Vaccination with 10-Valent Pneumococcal Non-Typeable Haemophilus influenza Protein D Conjugate Vaccine (PHiD-CV) on the Nasopharyngeal Microbiome of Kenyan Toddlers. PLOS ONE 10: e0128064. doi: 10.1371/journal.pone.0128064
[21]	Prevaes SM, de Winter-de Groot KM, Janssens HM, et al. (2015) Development of the Nasopharyngeal Microbiota in Infants with Cystic Fibrosis. Am J Respir Crit Care Med.
[22]	Cremers AJ, Zomer AL, Gritzfeld JF, et al. (2014) The adult nasopharyngeal microbiome as a determinant of pneumococcal acquisition. Microbiome 2: 44. doi: 10.1186/2049-2618-2-44
[23]	Allen EK, Koeppel AF, Hendley JO, et al. (2014) Characterization of the nasopharyngeal microbiota in health and during rhinovirus challenge. Microbiome 2: 22. doi: 10.1186/2049-2618-2-22
[24]	Biesbroek G, Bosch AA, Wang X, et al. (2014) The impact of breastfeeding on nasopharyngeal microbial communities in infants. Am J Respir Crit Care Med 190: 298–308.
[25]	Sakwinska O, Bastic Schmid V, Berger B, et al. (2014) Nasopharyngeal microbiota in healthy children and pneumonia patients. J Clin Microbiol 52: 1590–1594. doi: 10.1128/JCM.03280-13
[26]	Bassis CM, Tang AL, Young VB, et al. (2014) The nasal cavity microbiota of healthy adults. Microbiome 2: 27. doi: 10.1186/2049-2618-2-27
[27]	Perez-Losada M, Castro-Nallar E, Bendall ML, et al. (2015) Dual Transcriptomic Profiling of Host and Microbiota during Health and Disease in Pediatric Asthma. PLOS ONE 10: e0131819. doi: 10.1371/journal.pone.0131819
[28]	Castro-Nallar E, Shen Y, Freishtat RJ, et al. (2015) Integrating metagenomics and host gene expression to characterize asthma-associated microbial communities. BMC Medical Genomics 8: 50. doi: 10.1186/s12920-015-0121-1
[29]	Bogaert D, Keijser B, Huse S, et al. (2011) Variability and diversity of nasopharyngeal microbiota in children: a metagenomic analysis. PLOS ONE 6: e17035. doi: 10.1371/journal.pone.0017035
[30]	Pérez-Losada M, Crandall KA, Freishtat RJ (2016) Two sampling methods yield distinct microbial signatures in the nasopharynx of asthmatic children. Microbiome[in press].
[31]	Benton AS, Wang Z, Lerner J, et al. (2010) Overcoming heterogeneity in pediatric asthma: tobacco smoke and asthma characteristics within phenotypic clusters in an African American cohort. J Asthma 47: 728–734. doi: 10.3109/02770903.2010.491142
[32]	Kozich JJ, Westcott SL, Baxter NT, et al. (2013) Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl Environ Microbiol 79: 5112–5120. doi: 10.1128/AEM.01043-13
[33]	Schloss PD, Westcott SL, Ryabin T, et al. (2009) Introducing mothur: Open-Source, Platform-Independent, Community-Supported Software for Describing and Comparing Microbial Communities. Appl Environ Microbiol 75: 7537–7541. doi: 10.1128/AEM.01541-09
[34]	Schloss PD, Gevers D, Westcott SL (2011) Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies. PLOS ONE 6: e27310. doi: 10.1371/journal.pone.0027310
[35]	Edgar RC, Haas BJ, Clemente JC, et al. (2011) UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27: 2194–2200. doi: 10.1093/bioinformatics/btr381
[36]	Wang Q, Garrity GM, Tiedje JM, et al. (2007) Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol 73: 5261–5267. doi: 10.1128/AEM.00062-07
[37]	Caporaso JG, Kuczynski J, Stombaugh J, et al. (2010) QIIME allows analysis of high-throughput community sequencing data. Nat Methods 7: 335–336. doi: 10.1038/nmeth.f.303
[38]	Price MN, Dehal PS, Arkin AP (2010) FastTree 2-Approximately Maximum-Likelihood Trees for Large Alignments. PLOS ONE 5.
[39]	Faith DP (1992) Conservation evaluation and phylogenetic diversity. Biol Conserv 61: 1–10. doi: 10.1016/0006-3207(92)91201-3
[40]	Dixon P (2003) VEGAN, a package of R functions for community ecology. J Veg Sci 14: 927–930. doi: 10.1111/j.1654-1103.2003.tb02228.x
[41]	RStudioTeam (2015) RStudio: Integrated Development for R. RStudio, Inc., Boston, MA URL http://www.rstudio.com/.
[42]	Biesbroek G, Wang X, Keijser BJ, et al. (2014) Seven-valent pneumococcal conjugate vaccine and nasopharyngeal microbiota in healthy children. Emerg Infect Dis 20: 201–210. doi: 10.3201/eid2002.131220
[43]	Yan M, Pamp SJ, Fukuyama J, et al. (2013) Nasal microenvironments and interspecific interactions influence nasal microbiota complexity and S. aureus carriage. Cell Host Microbe 14: 631–640. doi: 10.1016/j.chom.2013.11.005
[44]	Cremers AJH, Zomer AL, Gritzfeld JF, et al. (2014) The adult nasopharyngeal microbiome as a determinant of pneumococcal acquisition. Microbiome 2.
[45]	Bassis CM, Erb-Downward JR, Dickson RP, et al. (2015) Analysis of the upper respiratory tract microbiotas as the source of the lung and gastric microbiotas in healthy individuals. MBio 6: e00037.
[46]	Morgan JL, Darling AE, Eisen JA (2010) Metagenomic sequencing of an in vitro-simulated microbial community. PLOS ONE 5: e10209. doi: 10.1371/journal.pone.0010209
[47]	Yuan S, Cohen DB, Ravel J, et al. (2012) Evaluation of methods for the extraction and purification of DNA from the human microbiome. PLOS ONE 7: e33865. doi: 10.1371/journal.pone.0033865
[48]	P OC, Aguirre de Carcer D, Jones M, et al. (2011) The effects from DNA extraction methods on the evaluation of microbial diversity associated with human colonic tissue. Microb Ecol 61: 353–362. doi: 10.1007/s00248-010-9771-x
[49]	Mackenzie BW, Waite DW, Taylor MW (2015) Evaluating variation in human gut microbiota profiles due to DNA extraction method and inter-subject differences. Front Microbiol 6: 130.

This article has been cited by:

1.	Laura Lomba, David Lapeña, Natalia Ros, Elena Aso, Mariachiara Cannavò, Diego Errazquin, Beatriz Giner, Ecotoxicological study of six drugs in Aliivibrio fischeri, Daphnia magna and Raphidocelis subcapitata, 2020, 27, 0944-1344, 9891, 10.1007/s11356-019-07592-8
2.	Laura Lomba, M. Pilar Ribate, Estefanía Zuriaga, Cristina B. García, Beatriz Giner, Acute and subacute effects of drugs in embryos of Danio rerio. QSAR grouping and modelling, 2019, 172, 01476513, 232, 10.1016/j.ecoenv.2019.01.081
3.	Carlos Castillo-Zacarías, Mario E. Barocio, Enrique Hidalgo-Vázquez, Juan Eduardo Sosa-Hernández, Lizeth Parra-Arroyo, Itzel Y. López-Pacheco, Damià Barceló, Hafiz N.M. Iqbal, Roberto Parra-Saldívar, Antidepressant drugs as emerging contaminants: Occurrence in urban and non-urban waters and analytical methods for their detection, 2021, 757, 00489697, 143722, 10.1016/j.scitotenv.2020.143722
4.	Melanie M. Marshall, Kevin E. McCluney, Mixtures of co-occurring chemicals in freshwater systems across the continental US, 2021, 268, 02697491, 115793, 10.1016/j.envpol.2020.115793
5.	Harveer S. Srain, Karen F. Beazley, Tony R. Walker, Pharmaceuticals and personal care products and their sublethal and lethal effects in aquatic organisms, 2021, 29, 1181-8700, 142, 10.1139/er-2020-0054
6.	Kaiwen Hu, Weifeng Li, Weixia Zhang, Kuankuan Yuan, Chenxin Gong, Yang Shu, Yingying Yu, Conghui Shan, Yan Gao, Xunyi Zhang, Haibin Yu, Wei Shi, Guangxu Liu, Diltiazem disrupts Ca2+-homeostasis and exerts immunotoxic effects on a marine bivalve mollusc, the blood clam (Tegillarca granosa), 2025, 217, 0025326X, 118055, 10.1016/j.marpolbul.2025.118055

Reader Comments

Your name:*

Email:*
© 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)