text:
From Genomics to Systems Biology
MICRAT: a novel algorithm for inferring gene regulatory networks using time series gene expression data
Bei Yang1,2,
Yaohui Xu1,
Andrew Maxwell3,
Wonryull Koh3,
Ping Gong4 &
Chaoyang Zhang3
BMC Systems Biology volume 12, Article number: 115 (2018) Cite this article
Reconstruction of gene regulatory networks (GRNs), also known as reverse engineering of GRNs, aims to infer the potential regulation relationships between genes. With the development of biotechnology, such as gene chip microarray and RNA-sequencing, the high-throughput data generated provide us with more opportunities to infer the gene-gene interaction relationships using gene expression data and hence understand the underlying mechanism of biological processes. Gene regulatory networks are known to exhibit a multiplicity of interaction mechanisms which include functional and non-functional, and linear and non-linear relationships. Meanwhile, the regulatory interactions between genes and gene products are not spontaneous since various processes involved in producing fully functional and measurable concentrations of transcriptional factors/proteins lead to a delay in gene regulation. Many different approaches for reconstructing GRNs have been proposed, but the existing GRN inference approaches such as probabilistic Boolean networks and dynamic Bayesian networks have various limitations and relatively low accuracy. Inferring GRNs from time series microarray data or RNA-sequencing data remains a very challenging inverse problem due to its nonlinearity, high dimensionality, sparse and noisy data, and significant computational cost, which motivates us to develop more effective inference methods.
We developed a novel algorithm, MICRAT (Maximal Information coefficient with Conditional Relative Average entropy and Time-series mutual information), for inferring GRNs from time series gene expression data. Maximal information coefficient (MIC) is an effective measure of dependence for two-variable relationships. It captures a wide range of associations, both functional and non-functional, and thus has good performance on measuring the dependence between two genes. Our approach mainly includes two procedures. Firstly, it employs maximal information coefficient for constructing an undirected graph to represent the underlying relationships between genes. Secondly, it directs the edges in the undirected graph for inferring regulators and their targets. In this procedure, the conditional relative average entropies of each pair of nodes (or genes) are employed to indicate the directions of edges. Since the time delay might exist in the expression of regulators and target genes, time series mutual information is combined to cooperatively direct the edges for inferring the potential regulators and their targets. We evaluated the performance of MICRAT by applying it to synthetic datasets as well as real gene expression data and compare with other GRN inference methods. We inferred five 10-gene and five 100-gene networks from the DREAM4 challenge that were generated using the gene expression simulator GeneNetWeaver (GNW). MICRAT was also used to reconstruct GRNs on real gene expression data including part of the DNA-damaged response pathway (SOS DNA repair network) and experimental dataset in E. Coli. The results showed that MICRAT significantly improved the inference accuracy, compared to other inference methods, such as TDBN, etc.
In this work, a novel algorithm, MICRAT, for inferring GRNs from time series gene expression data was proposed by taking into account dependence and time delay of expressions of a regulator and its target genes. This approach employed maximal information coefficients for reconstructing an undirected graph to represent the underlying relationships between genes. The edges were directed by combining conditional relative average entropy with time course mutual information of pairs of genes. The proposed algorithm was evaluated on the benchmark GRNs provided by the DREAM4 challenge and part of the real SOS DNA repair network in E. Coli. The experimental study showed that our approach was comparable to other methods on 10-gene datasets and outperformed other methods on 100-gene datasets in GRN inference from time series datasets.
Reconstruction of gene regulatory networks (GRNs), also known as reverse engineering [1] of GRNs, aims to infer the potential relationship between genes using gene expression data. Sometimes the expression of a gene is affected by others (named as its regulators) and meanwhile regulates its downstream target genes. Reconstruction of gene regulatory networks can help people understand the mechanism of gene interactions in the cell and reveal the mystery of life. The discovering of the regulation relationship between genes was initially achieved by biological experiments [2]. However, this approach was costly and progressed slowly in research, thus becoming a bottleneck that restricted the development of biological systems. With the development of biotechnology, such as gene chip microarray, the high-throughput data generated provides us with more opportunities to understand the potential regulatory relationships of genes. In recent years, developing machine learning and data mining algorithms and applying them to reconstruction of gene regulatory networks has become an active research topic in bioinformatics.
Different computational methods have been proposed to tackle the gene regulatory networks identification problems. The typical methods include regression models [3, 4], Bayesian networks [5,6,7,8,9,10], state space models [11,12,13,14] and information theory models [15,16,17,18,19,20], etc.
Several existing methods formulate gene regulatory networks as regression problems which are based on the assumption that all regulatory interactions are linear, meaning that the gene expression level of a target gene is a combination of the expression levels of its transcriptional factors. However, biological networks are known to exhibit multiple regulation mechanisms including non-linear interactions.
The Bayesian network model employed the joint probability distribution of gene expression data to construct a directed acyclic graph to represent the relationship between genes. A relative change ratio (RCR) method was proposed by Li et al. to preprocess the null-mutants steady state data in order to find the key genes and build GRNs, in which these selected key genes have a higher potential than other genes to play very critical roles [10]. The LBN algorithm (local Bayesian network) [21] was designed by Liu et al. to improve the accuracy of GRN inference from gene expression data by exploring advantages of Bayesian network (BN) and conditional mutual information (CMI) methods. The LBN algorithm first uses CMI to construct an initial network or GRN. Then, the BN method is employed to generate a series of local BNs by selecting the k-nearest neighbors of each gene as its candidate regulatory genes. Integrating these local BNs forms a tentative network by performing CMI. The final network, or GRN, can be obtained by iteratively performing CMI and local BN on the tentative network. TDBN [7] was used to infer GRNs from time series trajectory data which were combined with the previous knowledge gained in the above step. Wang et al. developed the pLasso [22] method from a Bayesian perspective for the reconstruction of gene networks. This method assigns different prior distributions to different edge subsets according to a modified Bayesian information criterion that incorporates prior knowledge on both the network structure and the pathway information. It used information from the Pathway Commons (PC) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) as prior information for the network reconstruction. However, in reality there is no prior knowledge while inferring the interaction relationship between genes. Moreover, signaling pathways are dynamic events that take place over a given period of time. In order to identify these pathways, expression data over time is required. Dynamic Bayesian network (DBN) is an important and widely used approach for predicting the GRNs from time series expression data but has high computational costs and may not be able to infer a large GRN.
Furthermore, other methods for inferring GRNs have been proposed. Patrilia et al. proposed a new algorithm, namely iRafNet [23], which integrated different data types under a unified random forest framework. The key idea of iRafNet is to introduce a weighted sampling scheme within random forest to incorporate information from other sources of data. Specifically, the model considers the expression of each gene as a function of the expression of other genes. In literature [24], a fast method was developed for inferring gene regulatory relationships from just knockdown data. It used a simple linear regression model focusing on single regulator-target gene pairs based on knockdown data, and allowed the incorporation of prior knowledge about the relationships and generates posterior probabilities.
Mutual information (MI) is another measurement of the dependence of two variables and is often used to assess the relationship between genes. Margolin [16] and Sales [17] employed MI between two genes to measure their associations in order to reconstruct gene regulatory networks. But they could not tackle the problem that one gene is regulated by multiple genes. PCA_CMI [18] was proposed by Zhang et al. for inferring GRNs from gene expression data considering the non-linear dependence and topological structure of GRNs by employing path consistency algorithm based on conditional mutual information (CMI). It discovered the associations between related genes without identifying the regulators and their targets.
Reshef et al. presented a measure of dependence for two-variable relationships, named the maximal information coefficient (MIC) [25], which captures a wide range of functional and nonfunctional associations. MIC has better performance on measuring the dependence of two-variable relationships in comparison with other methods. Thus, we concentrate our effort on identifying the interactions between genes based on this method.
In this paper, we developed a novel algorithm, MICRAT (Maximal Information coefficient with Conditional Relative Average entropy and Time series mutual information), for inferring GRNs from time series gene expression data. This approach employed MIC as measurement for assessing the correlation between each pair of genes while constructing an undirected graph to represent the underlying relationship between genes. It subsequently combined conditional relative average entropy with time series mutual information to determine the potential regulator and its target genes.
Information theory has been used to reconstruct GRNs from gene expression data. Mutual information is generally used as a useful criterion for measuring the dependence between variables (genes) X and Y [18]. For gene expression data, variable X is a vector, in which the element x denotes gene expression value in different conditions (samples/time points). In this section, we introduced several basic concepts of information theory and graph theory which will be used in the proposed method.
Mutual information and conditional mutual information
For two random variables (genes) X and Y, the mutual information of X and Y is defined as:
$$ I\left(X,Y\right)=H(X)+H(Y)-H\left(X,Y\right) $$
where H(X) and H(Y) are the entropy of the random variables X and Y, respectively; H(X, Y) is the joint entropy of X and Y.
For a discrete variable X, the entropy H(X) measures the average uncertainty of variable X, and can be computed by
$$ H(X)=-{\sum}_{x\in X}p(x)\log p(x) $$
where p(x) is the probability that the variable X takes x;
The joint entropy of X and Y is denoted by
$$ H\left(X,Y\right)=-{\sum}_{x\in X,y\in Y}p\left(x,y\right)\log p\left(x,y\right) $$
where p(x, y) is the joint probability of X and Y.
The mutual information of X and Y is calculated by
$$ I\left(X,Y\right)={\sum}_{x\in X,y\in Y}p\left(x,y\right)\log \frac{p\left(x,y\right)}{p(x)p(y)} $$
The conditional mutual information [16] of X and Y given Z, I(X, Y| Z) is defined as
$$ I\left(X,Y|Z\right)={\sum}_{x\in X,y\in Y,z\in Z}p\left(x,y,z\right)\log \frac{p\left(x,y|z\right)}{p\left(x|z\right)p\left(y|z\right)} $$
where p(x, y, z) is the joint probability of variables X, Y and Z; p(x| z) is the conditional probability that the variable x holds under the condition that the variable z is established; similarly, p(x, y| z) is the conditional probability that the variable x and the variable y are established under the condition that the variable z is satisfied.
In this paper, the probability of a variable (gene) is estimated with the Gaussian kernel probability density estimator [26] as follows [18].
$$ P\left({X}_i\right)=\frac{1}{n}{\sum}_{j=1}^n\frac{1}{{\left(2\pi \right)}^{m/2}{\left|C\right|}^{m/2}}\exp \left(-\frac{1}{2}{\left({X}_j-{X}_i\right)}^T{C}^{-1}\left({X}_j-{X}_i\right)\right) $$
where C is the covariance matrix of variable X, ∣C∣ is the determinant of matrix C, n is the number of samples (time points), and m is the number of variables (genes).
Then the entropy of variable X can be computed by using (2) and (6).
$$ H(X)=\log \left[{\left(2\pi e\right)}^{\frac{m}{2}}{\left|C\right|}^{\frac{1}{2}}\right]=\frac{1}{2}\log {\left(2\pi e\right)}^m\mid C\mid $$
With formulas (1) and (7), the mutual information of variables X and Y can be obtained as given in Eq. (8).
$$ I\left(X,Y\right)=\frac{1}{2}\log \frac{\mid C(X)\mid \bullet \mid C(Y)\mid }{\mid C\left(X,Y\right)\mid } $$
Similarly, the conditional mutual information of variables X and Y given variable Z can be calculated by Eq. (9).
$$ I\left(X,Y|Z\right)=\frac{1}{2}\log \frac{\mid C\left(X,Z\right)\mid \bullet \mid C\left(Y,Z\right)\mid }{\mid C(Z)\mid \bullet \mid C\left(X,Y,Z\right)\mid } $$
MI can be used for measuring the dependence of two variables (genes) and CMI can help determine the dependence between two variables given another variable. Generally speaking, a higher MI (or CMI) indicates closer relationship between variables.
Maximal information coefficient
The maximal information coefficient (MIC) is an effective measurement of interesting relationship between pairs of variables. MIC captures a wide range of functional and nonfunctional associations, and has two heuristic properties: generality and equitability. The generality ensures that with sufficient sample size MIC captures a wide range of interesting relationships, not limited to specific function types, or even to all functional relationships; and the equitability means that MIC gives similar scores to equally noisy relationships of different types. These two heuristic properties of MIC indicate that it is suitable to infer various underlying relationships.
The Maximal Information Coefficient of a set D of two-variable data with sample size n and grid size less than B(n) is given by literature [25].
$$ MIC(D)={\mathit{\max}}_{XY<B(n)}M{(D)}_{X,Y}\ {=\mathit{\max}}_{XY<B(n)}\frac{I^{\ast}\left(D,X,Y\right)}{\mathit{\log}\left(\mathit{\min}\left(X,Y\right)\right)} $$
where the dataset D is composed of the pairs of values <x, y>, and x ∈ X, y ∈ Y; M(D)X, Y is the characteristic matrix of dataset D of variable X and Y with entries
$$ M{(D)}_{x,y}=\frac{I^{\ast}\left(D,x,y\right)}{\log \left(\min \left(x,y\right)\right)} $$
where I∗(D, x, y) is the maximum mutual information achieved by x and y; n is size of D; log( min(x, y)) is used to normalize the maximum mutual information; B(n) is the grid size, and B(n) = n0.6, which was found to work well in practice [25].
Conditional relative average entropy
The conditional relative average entropy (CRAE) [27] is defined as:
$$ CRAE\left(X\to Y\right)=\frac{H\left(Y|X\right)}{H(Y)\bullet \mid Y\mid } $$
where ∣Y∣ is the number of values of variable Y, and H(Y) is the entropy of the variable Y; H(Y| X) is the conditional entropy of the variable Y under the condition of the given variable X.
Time series mutual information
In most cases, there is a time delay during the gene regulations or at least the gene expressions change simultaneously for regulatory and target genes. Thus the time series mutual information could contribute to the orientation of the regulation.
Suppose that we have T time points and the expression levels of N genes are measured at each time point. The time series gene expression data can be summarized as a T × N matrix X = (x1, x2,…, xT)T whose ith row vector xi = (xi1, xi2, …xiN)Tcorresponds to a gene expression level vector measured at timepoint i.
For given genes X and Y, the time series mutual information from X to Y during time points 1 to T is defined as
$$ {I}_T\left(X\to Y\right)=I\left({X}_{1\to \left(T- tdelay\right)},{Y}_{1+ tdelay\to T}\right) $$
where tdelay(0 ≤ tdelay ≤ T) represents the time delay for gene X regulating gene Y; X1 → (T − tdelay) is a vector representing the gene expression level of X from timepoint 1 to timepoint (T − tdelay). Similarly, Y1 + tdelay → T represents the gene expression level of Y from timepoint 1 + tdelay to timepoint T. In this paper, we set the time delay to 1 since from the observation of our experiments there is no significant improvement to the accuracy with different settings.
Gene regulatory networks
In this work, a gene regulatory network is described by a directed graph (digraph). A graph G is defined by the pair (V(G), E(G)), where V(G) denotes the set of vertices (nodes) and E(G) ⊆ V(G) × V(G) denotes the set of edges. In a digraph, an edge is defined by an ordered pair of vertices (x, y) denoting the edge direction, from vertex x pointing to vertex y. Here, the nodes represent genes while the edges describe the gene regulatory interactions. A directed edge (x, y) indicates that gene x regulates gene y, while an undirected edge (x, y) only indicates that there is an association between gene x and gene y.
Algorithm of MICRAT
We propose a novel algorithm called MICRAT for reconstructing gene regulatory networks using time series gene expression profiles. The algorithm mainly consists of two procedures: First, an undirected graph is generated representing the associations between genes based on the maximal information coefficient (MIC) of each pair of genes. Then, the conditional relative average entropy combining with time series mutual information is used to determine the directions of the edge in the undirected graph. The flowchart of the method for reconstructing GRN is shown in Fig. 1, and more details are described in the following subsections.
Reconstruction of GRN flowchart
Generating an undirected gene interaction graph
The types of relationship between genes vary diversely. It is difficult, even impossible to identify the gene associations with only one function. The maximal information coefficient could capture a wide range of dependence of two variables, both functional and nonfunctional. Thus it is suitable for measuring the associations between genes.
We refer to two gene expression data GX and GY as variables X and Y, respectively. The association between them is measured by
$$ {I}_D\left(X,Y\right)= MIC\left(X,Y\right) $$
where MIC(X, Y) denotes the maximal information coefficient of X and Y. A higher value of ID(X, Y) indicates a closer relationship between X and Y, vice versa. The edges in the undirected graph are added based on a given a threshold θ of ID(X, Y).
The generation of an undirected gene interaction graph includes the following five steps:
Initialize the gene regulatory network, GD, as an undirected graph with no vertices and no edges. There are n genes in the dataset D.
For a given time series gene expression dataset D, calculate ID(X, Y) for each pair of genes (x, y).
Rank the pairs of genes by descending order of ID(X, Y) obtained by step (2) to a list called RLI.
From the top of RLI, check each pair of genes (x, y), add an undirected edge between vertices x and y if ID(X, Y) ≥ θ where θ is a threshold and will be discussed in the experiment section.
Delete redundant edges. There might be some redundant edges in the undirected graph especially when there are triangle cycles. Conditional mutual information of genes is used to determine whether to delete edges or not. For example, there is a triangle cycle including three genes X, Y and Z in the graph since ID(X, Y), ID(X, Z) and ID(Y, Z) are not lower than θ and three edges between X, Y and Z are included. We calculated three conditional mutual information values of I(X, Y| Z), I(X, Z| Y) and I(Y, Z| X), and deleted edge (Y, Z) if I(Y, Z| X) = 0, which indicates that gene Y has no relationship with gene Z.
Data discretization and edge-orientation
In our algorithm, the edges in the undirected graph are oriented by combining the conditional relative average entropy with time series mutual information of each pair of nodes. From the observation of our experiments, better performance was obtained on discretized data. Therefore, the gene expression data were discretized before being used for edge orientation.
The gene expression data are first normalized by standard fraction Z-Score and then discretized according to a given threshold. The standard fraction Z _ score of the gene X at the time point tj is defined as:
$$ {Z}_{x,j}=\frac{\mid {\mathrm{x}}_j-\mu \mid }{\sigma } $$
where μ denotes the mean value of the gene expression data at all time points of gene X, σ is the standard deviation, and xj denotes the expression value of the gene X at the time point tj. Given a threshold k, if Zx, j ≥ k, then the gene expression data at time tj is denoted as 1; otherwise the gene expression data is denoted as 0.
All edges in the undirected graph are oriented by the following procedure: Given two vertices X and Y of an edge in the undirected graph, if CRAE(X → Y) + IT(X → Y) > CRAE(Y → X) + IT(Y → X), then the orientation is X → Y, indicating that gene X regulates gene Y; While if CRAE(X → Y) + IT(X → Y) < CRAE(Y → X) + IT(Y → X), then the orientation is Y → X, which means gene Y regulates gene X. In the case of CRAE(X → Y) + IT(X → Y) = CRAE(Y → X) + IT(Y → X), the higher value of IT(X → Y) (or IT(Y → X)) indicates X → Y (or Y → X). The pseudo code of algorithm MICRAT is given in Fig. 2.
Pseudo code of algorithm MICRAT
In this section, we evaluate the performance of algorithm MICRAT by applying it to synthetic datasets as well as a real gene expression dataset for inferring GRNs. Several experiments were carried out to demonstrate the performance of our approach by comparison with other related methods.
For synthetic datasets, five 10-genes and five 100-genes time series expression datasets from the DREAM4 challenge were used, respectively [28]. The real life gene expression dataset involves an eight-gene subnetwork, part of a DNA-damage response pathway (SOS pathway) in the bacteria E. Coli.
Generally, there are several measurements for evaluating the performance of GRN inference methods, including Precision, Recall, Accuracy (ACC), F _ score, Matthews correlation coefficient (MCC), etc. They are defined as follows:
$$ Precision=\frac{TP}{TP+ FP} $$
$$ Recall=\frac{TP}{TP+ FN} $$
$$ ACC=\frac{TP+ TN}{TP+ FP+ TN+ FN} $$
$$ F\_ score=\frac{2\times Precision\times Recall}{Precision+ Recall} $$
$$ MCC=\frac{TP\times TN- FP\times FN}{\sqrt{\left( TP+ FP\right)\left( TP+ FN\right)\left( TN+ FP\right)\left( TN+ FN\right)}} $$
where TP, FP, TN, FN are the numbers of true positives, false positives, true negatives and false negatives, respectively. In this study, true positives indicate correctly inferred edges which exist in golden standard networks, while false positives denote wrongly inferred edges which do not exist in golden standard networks, and so on for true negatives and false negatives, respectively. Usually, there exists an inverse relationship between precision and recall. It is possible to increase one at the cost of reducing the other. Therefore, Accuracy, F _ score and Matthews correlation coefficient are widely used as balanced evaluation measures. Here we use ACC, F _ score and MCC as measurements for evaluating the performance of our method.
In order to set the best value of θ, the threshold of MIC of two gene expression data while measuring the independence between genes, we test MICRAT by varying θ from 0 to 1 with interval of 0.01 on both 10-gene expression datasets and 100-gene expression datasets. Empirical results showed that with the increase of θ, the precision of the model becomes higher while the recall goes down. According to the experimental results, the comparatively good performance of MICRAT could be obtained while the thresholds were set to 0.43 for 10-gene datasets and 0.15 for 100-gene datasets, respectively.
Experiments on simulation datasets
We first evaluated MICRAT on simulation datasets which were generated based on benchmarking networks. Many tools have been developed for assessing the effectiveness of GRN inference methods. The DREAM4 challenge introduced a critical performance assessment framework of methods for GRN inference and presented an in silico benchmark suite as a blinded, community-wide challenge within the context of the DREAM project. In this challenge, the gene expression datasets with noise and their golden standard (benchmark) networks were given. The golden standard networks were selected from source networks of real species [28]. They are widely used as a benchmark for the evaluation of GRN inference methods. Here, we tested our algorithm on the DREAM4 challenge time series gene expression data in networks of size 10 and 100 genes, respectively.
As mentioned above, MICRAT found the associations between pairs of genes in the beginning. Then the conditional relative average entropy in cooperation with time series mutual information was used to orient the edges in the undirected graph and obtain the GRN. For evaluating the performance of our approach, we compared our algorithm with TDBN [7, 10] and NARROMI [19]. TDBN is a typical method with significant impact on inferring GRNs from time series gene expression data. NARROMI is quite related to our method and it significantly outperformed other methods, such as LP [29], LASSO [30, 31], ARACNE [16] and GENIE3 [32]. This procedure was tested on DREAM4 time series gene expression datasets with 10-gene networks and 100-gene networks, respectively. Experimental results were given in Table 1. First, we carried out MICRAT on five 10-gene datasets. There are 21 × 5 = 105 time points in each of the datasets. We computed ACC, F _ score and MCC on each of the five datasets as evaluation measures for the compared methods. Here, we chose 0.43 as the threshold value of MIC to determine the dependence between genes in the inferred network for our method, and the threshold of Z _ score was set to k = 1.2 for discretization of gene expression data when the edge is oriented. For NARROMI, the threshold is set to 0.05 which is the default value given in literature [19]. For TDBN, all the parameters were chosen as default in their literature. The accuracy of our approach is higher than those of the other methods in four datasets and slightly lower than that of NARROMI for dataset 5. The left part of Fig. 3 shows the average accuracies on five 10-gene datasets for three methods. It can be observed that MICRAT performed better than the others. For the F _ score measurement, MICRAT has the best performance on two datasets, but the worst performance on one dataset. Figure 4 shows that the average F _ scores on five 10-gene datasets is somewhat inferior to NARROMI but superior to TDBN. From Table 1 it can be observed that MICRAT is only slightly inferior to NARROMI on one dataset but has better performance on others for MCC. Figure 5 gives the average Matthews correlation coefficients of the three methods on five 10-gene datasets and five 100-gene datasets, respectively. We can see that MICRAT is superior to NARROMI and significantly outperforms TDBN.
Table 1 Experimental results for three methods on both 10-gene datasets and 100-gene datasets
Comparison of three methods for average ACC on 10-gene and 100-gene datasets
Comparison of three methods for average F_score on 10-gene and 100-gene datasets
Comparison of three methods for average MCC on 10-gene and 100-gene datasets
We also applied our algorithm to the five DREAM4 100-gene datasets. Each dataset has 21 × 10 = 210 time points. We also computed ACC, F _ score and MCC on each of the five datasets as evaluation measures for MICRAT and the other two methods. Experimental results were showed in Table 1 and the right part of Figs. 3, 4 and 5. In these experiments, we chose 0.15 as the threshold of MIC to determine the regulation relationship between genes in the inferred network. For NARROMI, the threshold is also set to 0.05, the same value as in literature [19]. As shown in Table 1 and the right part of Figs. 3, 4 and 5, all the three measures of ACC, F _ score and MCC in our method are superior to NARROMI and TDBN on each of the five 100-gene time series dataset. The results show that our approach has better performance than other methods on the 100-gene datasets. Thus MICRAT has the potential to infer large GRNs using time series expression data with a large number of genes.
Experiments on real gene expression dataset
To validate our method on a real biological gene regulatory network, we analyzed the well-known SOS DNA repair network in E. Coli. This GRN is well known for its responsibility of repairing the DNA if it gets damaged. It is the largest, most complex and best understood DNA damage-inducible network to be characterized to date [33]. We applied our MICRAT algorithm to an eight-gene network, part of the SOS DNA repair network in the bacteria E. coli [34,35,36]. The expression data sets of the network were obtained from Uri Aron Lab [36]. These data are kinetics of 8 genes namely uvrD, lexA, umuD, recA, uvrA, uvrY, ruvA and polB. Four experiments were done at various UV light intensities (Exp.1 and 2: 5Jm− 2, Exp. 3 and 4: 20Jm− 2). In each experiment, the 8 genes were monitored at 50 instants which are evenly spaced by 6 min intervals. In order to assess the effectiveness of our algorithm, we compared the inferred network with the known interactions between these eight genes. In our experiment, the threshold of MIC was set to θ=0.1 for measuring the dependence of the pair of genes, and the threshold of Z _ score was also set to k = 1.2 for discretization of gene expression data while the edge is oriented. Figure 6 gives the real SOS DNA repair network. Figure 7 shows part of inferred GRN inferred by our approach. In this Figure, one can see that our method finds 6 out of the 9 edges in the target network and identifies lexA as the 'hub' gene for this network. The exact ground truth for this network is not precisely known, hence it is not possible to calculate well-known performance measures [33]. By using the known interactions obtained from literature [33, 37, 38], Table 2 showed the comparison of our algorithm with other methods such as Perrin [39], BANJO [40], GlobalMIT [41] and Morshed [33] for correct predictions on the regulatory relationships between these eight genes. From Table 2, it could be observed that our method identified most interactions except lexA→ruvA, lexA→lexA and recA→lexA, where the latter two are not inferred since our method does not consider the double regulation between two genes and the gene that is self-regulated. That is the limitation of our method, which motivates us to continue further study for improvement.
The target SOS DNA repair network
Reconstruction of SOS DNA repair network by MICRAT
Table 2 Comparison of interactions inferred by MICRAT with other methods on SOS DNA dataset
In this work, a novel algorithm, MICRAT, for inferring GRNs from time series gene expression data was proposed by taking into account dependence and time delay of expressions of a regulator and its target genes. This approach employed maximal information coefficients for reconstructing an undirected graph to represent the underlying relationships between genes. The edges were directed by combining conditional relative average entropy with time course mutual information of pairs of genes. The performance of our proposed algorithm was evaluated using synthetic datasets from benchmark GRNs provided by the DREAM4 challenge and the real, experimental dataset SOS DNA repair network in E. Coli. Experimental study showed that our approach was comparable to other methods on 10-gene datasets and outperformed other methods on 100-gene datasets in GRN inference from time series datasets. In the follow-up study, we plan to test the algorithm on large-scale datasets, including time course mRNA-Seq, human gene expression data, etc. and compare our approach with more recently proposed algorithms so as to evaluate its robustness and improve its performance.
ACC:
CMI:
Conditional mutual information
CRAE:
GRN:
MCC:
Matthews correlation coefficient
MIC:
Vijesh N, Chakrabarti SK, Sreekumar J. Modeling of gene regulatory networks: a review. J Biomed Sci Eng. 2013;06(1):223–31. http://file.scirp.org/pdf/JBiSE_2013022716483315.pdf. Accessed 03 July 2018.
Lemmens K, Dhollander T, Bie TD, et al. Inferring transcriptional modules from ChIP-chip, motif and microarray data. Genome Biol. 2006;7(5):R37. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.279.1345&rep=rep1&type=pdf. Accessed 03 July 2018.
PubMed PubMed Central Article CAS Google Scholar
Singh N, Vidyasagar M. bLARS: an algorithm to infer gene regulatory networks. IEEE/ACM Trans Comput Biol Bioinform. 2016;13(2):1. https://www.computer.org/csdl/trans/tb/2016/02/07138615.pdf. Accessed 03 July 2018.
Haury AC, et al. TIGRESS: trustful inference of gene regulation using stability selection. BMC Syst Biol. 2012;6(1):145. https://core.ac.uk/download/pdf/51228782.pdf. Accessed 03 July 2018.
PubMed PubMed Central Article Google Scholar
Friedman N, Linial M, Nachman I, et al. Using Bayesian network to analyze expression data. J Comput Biol. 2000;7(3–4):601–20. www.cs.huji.ac.il/~nir/Papers/FLNP1Full.pdf. Accessed 03 July 2018.
Murphy K, Mian S. Modeling gene expression data using dynamic Bayesian networks. Technical Report, Computer Science Division. Berkeley: University of California; 1999. www.cs.ubc.ca/~murphyk/Papers/ismb99.pdf. Accessed 03 July 2018
Zou M, Conzen SD. A new dynamic Bayesian network (DBN) approach for identifying gene regulatory networks from time course microarray data. Bioinformatics. 2005;21(1):71–9. https://academic.oup.com/bioinformatics/article/21/1/71/212416. Accessed 03 July 2018.
CAS PubMed Article Google Scholar
Yeung KY, Dombek KM, Lo K, et al. Construction of regulatory networks using expression time-series data of a genotyped population. Proc Natl Acad Sci U S A. 2011;108(48):19436–41. www.pnas.org/content/108/48/19436.full. Accessed 03 July 2018.
CAS PubMed PubMed Central Article Google Scholar
Haoni L, et al. Learning the structure of gene regulatory networks from time series gene expression data. BMC Genomics. 2011;12(Suppl 5):S13. https://www.ncbi.nlm.nih.gov/pubmed/22369588. Accessed 03 July 2018.
Li P, et al. Gene regulatory network inference and validation using relative chang tratio analysis and time-delayed dynamic Bayesian network. EURASIP J Bioinform Syst Biol. 2014;2014:12. https://link.springer.com/article/10.1186/s13637-014-0012-3. Accessed 03 July 2018.
Wu X, et al. State space model with hidden variables for reconstruction of gene regulatory networks. BMC Syst Biol. 2011;5(Suppl 3):S3. http://europepmc.org/articles/PMC3287571. Accessed 03 July 2018.
Wu FX. Gene regulatory network modelling: a state-space approach. Int J Data Min Bioinform. 2008;2(1):1–14. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.192.4775&rep=rep1&type=pdf. Accessed 03 July 2018.
PubMed Article Google Scholar
Osamu H, Ryo Y, Seiya I, Rui Y, Tomoyuki H, Charnock-Jones DS, Cristin P, Satoru M. Statistical inference of transcriptional module-based gene networks from time course gene expression profiles by using state space models. Bioinformatics. 2008;24:932–42. https://academic.oup.com/bioinformatics/article/24/7/932/295736. Accessed 03 July 2018.
Article CAS Google Scholar
Kojima K, Rui Y, Seiya I, Mai Y, Masao N, Ryo Y, Teppei S, Kazuko U, Tomoyuki H, Noriko G, Satoru M. A state space representation of VAR models with sparse learning for dynamic gene networks. Genome Inform. 2009;22:56–68. https://www.jsbi.org/pdfs/journal1/IBSB09/IBSB09006.pdf. Accessed 03 July 2018.
Wang J, Chen B, Wang Y, et al. Reconstructing regulatory networks from the dynamic plasticity of gene expression by mutual information. Nucleic Acids Res. 2013;41(8):395–408. https://academic.oup.com/nar/article/41/8/e97/2409387. Accessed 03 July 2018.
Margolin AA, Nemenman I, Basso K, et al. ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics. 2006;7(1):S7. https://arxiv.org/pdf/q-bio/0410037. Accessed 03 July 2018.
Sales G, Romualdi C. Parmigene--a parallel R package for mutual information estimation and gene network reconstruction. Bioinformatics. 2011;27(13):1876–7 (2). https://academic.oup.com/bioinformatics/article/27/13/1876/184634. Accessed 03 July 2018.
Zhang X, Zhao XM, He K, et al. Inferring gene regulatory networks from gene expression data by path consistency algorithm based on conditional mutual information. Bioinformatics. 2012;28(1):98–104. https://dl.acm.org/citationcfm?id=2139373 Accessed 03 July 2018.
Zhang X, et al. NARROMI: a noise and redundancy reduction technique improves accuracy of gene regulatory network inference. Bioinformatics. 2013;29(1):106–13. https://core.ac.uk/download/pdf/52414020.pdf. Accessed 03 July 2018.
Chaitankar V, Ghosh P, Perkins EJ, et al. A novel gene network inference algorithm using predictive minimum description length approach. BMC Syst Biol. 2010;4(Suppl 1):S7. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2880413. Accessed 03 July 2018.
Liu F, Zhang SW, Guo WF, Wei ZG, Chen L. Inference of gene regulatory network based on local Bayesian networks. PLoS Comput Biol. 2016;12(8):e1005024. https://doi.org/10.1371/journal.pcbi.1005024 http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005024. Accessed 03 July 2018.
Wang Z. Incorporating prior knowledge into gene network study. Bioinformatics. 2013;29(20):2633–40. https://academic.oup.com/bioinformatics/article/29/20/2633/277562. Accessed 03 July 2018.
Petralia F, Wang P, Yang J, Zhidong T. Integrative random forest for gene regulatory network inference. Bioinformatics. 2015;31:197–205. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4542785/. Accessed 03 July 2018.
Young WC, Raftery AE, Yeung KY. A posterior probability approach for gene regulatory network inference in genetic perturbation data. Math Biosci Eng. 2017;13(6):1241–51. https://arxiv.org/abs/1603.04835. Accessed 03 July 2018.
Reshef DN, Reshef YA, Finucane HK, et al. Detecting novel associations in large data sets. Science. 2011;334(6062):1518–24. http://science.sciencemag.org/content/334/6062/1518. Accessed 03 July 2018.
Basso K, et al. Reverse engineering of regulatory networks in human B cells. Nat Genet. 2005;37:382–92. https://s3-us-west-2.amazonaws.com/oww-files-public/b/b2/Basso.pdf. Accessed 03 July 2018.
Jiang J, Wang J, Yu H, et al. Poison identification based on Bayesian network: a novel improvement on K2 algorithm via Markov blanket[M]// advances in swarm intelligence: Springer Berlin Heidelberg; 2013. p. 173–82. https://link.springer.com/chapter/10.1007/978-3-642-38715-9_21. Accessed 03 July 2018
Dream4 In Silico Network Challenage. http://dreamchallenges.org/project/dream4-in-silico-network-challenge/. Accessed 03 July 2018.
Wang Y, Joshi T, Zhang XS, et al. Inferring gene regulatory networks from multiple microarray datasets. Bioinformatics. 2006;22(22):2413–20. https://academic.oup.com/bioinformatics/article/22/19/2413/240982. Accessed 03 July 2018.
Tibshirani R. Regression shrinkage and selection via the lasso. J Royal Stat Soc. 1996;58(1):267–88. http://statweb.stanford.edu/~tibs/ftp/lasso-retro.pdf. Accessed 03 July 2018.
Geeven G, et al. Identification of context-specific gene regulatory networks with GEMULA-gene expression modeling using LAsso. Bioinformatics. 2012;28(2):214–21. https://academic.oup.com/bioinformatics/article/28/2/214/198473. Accessed 03 July 2018.
Huynh-Thu VA, et al. Infering regulatory networks from expression data using tree-based methods. PLoS One. 2010;5(9):e12776. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0012776. Accessed 03 July 2018.
Morshed, et al. Simultaneous learning instantaneous and time-delayed genetic interactions using novel information theoretic scoring technique. BMC Syst Biol. 2012;6:62. https://link.springer.com/chapter/10.1007%2F978-3-642-24958-7_29. Accessed 03 July 2018.
Ronen M, Rosenberg R, Shraiman BI, et al. Assigning numbers to the arrows: parameterizing a gene regulation network by using accurate expression kinetics. Proc Natl Acad Sci U S A. 2002;99(16):10555–60. www.pnas.org/lookup/doi/10.1073/pnas.152046799. Accessed 03 July 2018.
Shenorr SS, Milo R, Mangan S, et al. Network motifs in the transcriptional regulation network of Escherichia coli. Nat Genet. 2002;31(1):64–8. https://www.weizmann.ac.il/mcb/UriAlon/sites/mcb.UriAlon/files/network_motifs_in_coli_0.pdf. Accessed 03 July 2018.
Uri Alon's SOS Dataset webpage. https://www.weizmann.ac.il/mcb/UriAlon/download/downloadable-data. Accessed 03 July 2018.
Noman N, Iba H. Inferring gene regulatory networks using differential evolution with local search heuristics. IEEE/ACM Trans Comput Biol Bioinform. 2007;4(4):643–7. https://www.ncbi.nlm.nih.gov/pubmed/17975274. Accessed 03 July 2018.
Kamura S, et al. Inference of S-system models of genetic networks using a cooperative coevolutionary algorithm. Bioinformatics. 2005;21(7):1154–63. https://academic.oup.com/bioinformatics/article/21/7/1154/268773. Accessed 03 July 2018.
Perrin BE, Ralaivola L, Mazurie A, et al. Gene networks inference using dynamic Bayesian networks. Bioinformatics. 2003;19(suppl 2):ii138–48. https://academic.oup.com/bioinformatics/article/19/suppl_2/ii138/180436. Accessed 03 July 2018.
Yu J, Smith VA, Wang PP, et al. Advances to Bayesian network inference for generating causal networks from observational biological data. Bioinformatics. 2004;20(18):3594–603. https://www.semanticscholar.org/paper/Advances-to-Bayesian-network-inference-for-causal-Yu-Smith/adcbcb725490d64d12b4f795e1e381ca6b8de4b4. Accessed 03 July 2018.
Vinh NX, Chetty M, Coppel R, et al. GlobalMIT: learning globally optimal dynamic bayesian network with the mutual information test criterion. Bioinformatics. 2011;27(19):2765–6. https://academic.oup.com/bioinformatics/article/27/19/2765/231220. Accessed 03 July 2018.
This work was partially supported by a subaward of the National Science Foundation EPS 0903787 and the U.S. Army Environmental Quality/Installation (EQI) Basic Research Program. Permission was granted by the Chief of Engineers to publish this information. It was also partially supported by the Key Research Projects of Universities in Henan Province of China under grant number 17A520015. The publication cost of this article was funded by the U.S. Army Environmental Quality/Installation (EQI) Basic Research Program.
The datasets used and/or analysed during the current study are available at http://dreamchallenges.org/project/dream4-in-silico-network-challenge/. The codes of the algorithm of MICRAT are available from the corresponding author on reasonable request.
About this supplement
This article has been published as part of BMC Systems Biology Volume 12 Supplement 7, 2018: From Genomics to Systems Biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume-12-supplement-7.
School of Information & Engineering, Zhengzhou University, Zhengzhou, 450000, China
Bei Yang & Yaohui Xu
Center of Precision Medicine, Zhengzhou University, Zhengzhou, 450000, China
Bei Yang
School of Computing, University of Southern Mississippi, Hattiesburg, MS, 39406, USA
Andrew Maxwell, Wonryull Koh & Chaoyang Zhang
Environmental Lab, US Army Engineer Research and Development Center, Vicksburg, MS, 39180, USA
Ping Gong
Yaohui Xu
Wonryull Koh
Chaoyang Zhang
BY, CZ and PG conceived the project. YX and BY developed and implemented the MICRAT algorithms and conducted the experimental study. BY and CZ wrote the paper. PG, AM and WK participated in the analysis and provided suggestions for the document. All authors read and approved the final manuscript.
Correspondence to Bei Yang or Chaoyang Zhang.
Yang, B., Xu, Y., Maxwell, A. et al. MICRAT: a novel algorithm for inferring gene regulatory networks using time series gene expression data. BMC Syst Biol 12, 115 (2018). https://doi.org/10.1186/s12918-018-0635-1
Time-series mutual information


All Contractors used by Bath and North East Somerset Council must first qualify to be on our framework of preferred suppliers.
If you are a contractor and wish to be considered, you must first register at ww.servingthesouthwest.org.uk. When Bath and North East Somerset are looking for new contractors, the opportunity will be published on this portal. If you are successful on qualifying for the framework, you will then be notified of all routes that are being tendered. For further information on the expectations and terms and conditions please see the Instructions to Tenderers and Provider's Handbook.
All drivers on B&NES Council contracts are not only representing themselves, but their company and Bath and North East Somerset Council. We have high expectations of all of our drivers and any allegation against them will be taken very seriously and investigated. As a guide for drivers, a handbook has been published with straightforward do's and don'ts and some helpful information on who to contact in emergencies and hidden disabilities. All companies awarded transport contracts by the Council are responsible for ensuring that the drivers they utilise on their transport have read and understood this document.
How long does a Driver's Licence last for?


In the nearly one year that the AMDVLK official Vulkan driver has been open-source there hasn't been any "releases" but rather new code drops on a weekly basis that is pushed out of their internal development repositories. But surprisingly this morning is now a v2018.4.1 release tag for this open-source AMD Vulkan Linux driver.
The AMDVLK public source repositories have just been perpetual Git while AMD pulls from their internal repositories when building out their official closed-source Windows/Linux Radeon Software driver releases (that also use their closed-source shader compiler currently rather than the open-source AMDGPU LLVM back-end, as used by the public AMDVLK sources). Waking up this morning there is now the first release tag in AMDVLK as v2018.4.1.
From the release is now a blessed source code drop as well as an Ubuntu 18.04 binary. There doesn't appear to be any real changes over the AMDVLK code drop from last week.
Perhaps this new release tag is to align with the Q4 Radeon Software driver releases -- which would be nice to know for comparison purposes as up to now that hasn't really been communicated. We should hear more soon and of course I'll work on firing up some benchmarks of it against the more common Mesa RADV Radeon Vulkan driver. This also marks the first official binary builds we are seeing of AMDVLK albeit only for Ubuntu 18.04 as generally users wanting to try out this driver have had to build it from source unless on Arch or the few other distributions where there are third-party package builds of it.
It will be interesting to see how AMDVLK evolves moving into 2019 considering that most Linux gamers remain devoted RADV users thanks to it generally working out well, being part of Mesa and so it's packaged much more widely, and generally "just works" well for Linux gaming.


This defaults to the top bank account in your chart of accounts, but you can always change it in the dropdown list to the right of the Deposit To box You can change the default account for Make Deposits, as well as Write Checks, Pay Bills and Pay Sales Tax by at Edit > Preferences > Checking.... Learn about the required links needed to be set up with Quickbooks for payroll to post properly. Changing the Bank Account Payroll Checks are posted to Learn how to change bank account payroll checks to post to a new one.
Deposit to: Choose the USD bank account Wells Fargo Bank USD Amount received: Matching the bank amount that was received Note: This will match in your online bank feeds if you have your bank account linked to QuickBooks Online.... This defaults to the top bank account in your chart of accounts, but you can always change it in the dropdown list to the right of the Deposit To box You can change the default account for Make Deposits, as well as Write Checks, Pay Bills and Pay Sales Tax by at Edit > Preferences > Checking.
Learn about the required links needed to be set up with Quickbooks for payroll to post properly. Changing the Bank Account Payroll Checks are posted to Learn how to change bank account payroll checks to post to a new one.... Click on the bank account/credit card and rename both the account code and the account name by adding an�A� in the end of the lines. (Xero doesn�t allow you to use the same account code and name on two accounts so therefore we need to rename it �A� =Archive) Press �Save�.
This defaults to the top bank account in your chart of accounts, but you can always change it in the dropdown list to the right of the Deposit To box You can change the default account for Make Deposits, as well as Write Checks, Pay Bills and Pay Sales Tax by at Edit > Preferences > Checking.
Click on the bank account/credit card and rename both the account code and the account name by adding an�A� in the end of the lines. (Xero doesn�t allow you to use the same account code and name on two accounts so therefore we need to rename it �A� =Archive) Press �Save�.


text:
Check out these 10 jobs of the future
by Gordon Hunt
9 Aug 20162.99k Views
A new report has found that virtual habitat designers, ethical technology advocates and freelance biohackers will be three of the jobs of the future.
When the Berlin Wall came down 27 years ago, on 9 November 1989, there was no such thing as a web developer.
When Michelle Smith won her third gold medal 20 years ago at the Atlanta 1996 Olympic Games, big data analysts were absent from the workforce.
Only a decade ago, in 2006, as Ireland played its final competitive rugby match at the old Lansdowne Road, social media managers were yet to be born.
Jobs are ever-changing, never more so than in the fast-paced digital revolution of recent decades. So Microsoft's paper on 10 future jobs makes for interesting reading.
Claiming that two-thirds of graduates in 2025 will be doing jobs that don't even exist yet, the Tomorrow's Jobs report – which The Future Laboratory and Microsoft collaborated on – predicts humanity will be curing disease and building virtual reality worlds from "the comfort of our own sofas".
Its 10 future jobs (the latter half of which are unlikely to be a reality in 2025), including some helpful descriptions, are as follows:
1. Virtual reality habitat designer
"By 2025, virtual habitat design will offer some of the most exciting and creative career prospects in a global industry that will be producing millions of new jobs.
"VR habitat designers will need to possess the storytelling skills of an online game designer and editor together with the spatial design expertise of an architect or town planner to be able to imagine and create entire virtual worlds."
2. Ethical technology advocate
"Ethical technology advocates will be mankind's go-betweens, with a wave of robots and artificial intelligence applications that will be helping to run our complex and connected world by 2025.
"One of their key jobs will be to negotiate our delicate relationship with the robots by setting the moral and ethical rules under which the machines – and their makers – operate and exist."
3. Digital cultural commentator
"In the 2020s, digital cultural commentator will be the secret weapon that both brands and centres of high culture will use to cut through the cacophony of online white noise to talk effectively to tomorrow's audience.
"Masters of the next generation of visual social media, they will be able to bridge the gap between the arts and a digitally savvy public by using simple, impactful images to communicate complex and challenging ideas."
Freelancer, via Shutterstock
4. Freelance biohacker
"Working from home, or from the growing number of freelance work hubs, freelance biohackers will work on open-source software platforms with hundreds, even thousands, of others in hive-like teams.
"University research departments and major drug and bioscience companies will use them to piece together complex DNA-based answers to some of the big questions of the next decade, from treatments for cancers in ageing populations to vaccines for new epidemics fuelled by our globalised culture and accelerating climate change."
5. Internet-of-things data creative
"It will be their job to sift through the waves of data being generated each day by devices in our clothes, our homes, our cars and our offices and find meaningful and useful ways to tell us what all that information is saying.
"They will need to have three key talents: a finely-honed ability to recognise patterns, a skill at asking sharp and difficult questions, and a natural flair for storytelling."
Space tour guide, via Shutterstock
6. Space tour guide
"Space tour guides will use their exhaustive knowledge of the location of the thousands of satellites and pieces of junk from previous space missions to construct visits to the most interesting places in orbit.
"They will have a secondary role in finding and plotting the location of forgotten abandoned spacecraft and defunct equipment to create a log of the orbiting memorabilia of mankind's first half century outside of Earth's atmosphere."
7. Personal content creator
"By the late 2020s, software-brain interfaces, pioneered by teams of neuroscientists, will have started to enter the mainstream, allowing mass audiences to read and capture thoughts, memories and dreams.
"Personal content creators will help people to use these systems to increase the storage capacity of their overstretched minds, providing services that allow them to dip in and out of treasured memories and experiences at will."
8. Rewilding strategist
"Rewilding strategists will stitch together viable ecosystems in stressed landscapes, using patchworks of flora and fauna from all over the world, rather than worrying about only using indigenous species.
"They will reintroduce plants and animals that have been extinct in a region for centuries (think wolves in Ireland) and manage assisted migrations in order to create resilient and vibrant landscapes in the face of advancing climate change."
9. Sustainable power innovator
"Sustainable power innovator – workers with expertise in chemistry and material science alongside finely-honed entrepreneurial instincts – will scour the periodic table, combining elements and organic materials to invent new battery storage capabilities, for both on and off-grid use.
"They will also oversee the introduction of super-fast charging facilities to cope with the power demands of our ever-growing reliance on the internet of things in a hyper-urbanised world."
Human body designer, via Shutterstock
10. Human body designer
"Over the next two decades, bio-engineering advances will extend the average healthy human lifespan to more than 100 years as the growth of replacement tissues and organs becomes every day and affordable.
"Human body designers will combine design skills with bio-engineering know-how to create a huge range of customised human limbs, either to perfectly match the existing skin tone, musculature and colour of the rest of a person's body, to provide exotic new looks or enhanced functionality for particular jobs or sports."
Looking for jobs in tech or science? Check out our Employer Profiles for information on companies hiring right now.
Main virtual reality image via Shutterstock
Related: future, biotech, Microsoft, internet of things, jobs, space, VR
Gordon joined Silicon Republic in October 2014 as a journalist, moving on to a new position as senior communications and content executive at NDRC in August 2017. Unafraid of heights or spiders, Gordon spends most of his time avoiding conversations about music, appreciating even the least creative pun and rueing the day he panicked when meeting Paul McGrath. His favourite thing on the internet remains the 'Random Article' link on Wikipedia.
'We've transformed our recruitment to ensure it reflects our values'
3 ways extended reality can help your employees

redpajama_set_name: RedPajamaCommonCrawl


text:
require_relative '../../spec_helper'
require_relative 'shared/new'

describe "IO.for_fd" do
  it_behaves_like :io_new, :for_fd
end

describe "IO.for_fd" do
  it_behaves_like :io_new_errors, :for_fd
end

redpajama_set_name: RedPajamaGithub


text:
\section{Introduction} \label{intro}

Ultracold and quantum degenerate molecular gases are currently a focus of interest within atomic physics, molecular physics and physical chemistry communities~\cite{Ni2010}.  The enormous control and high tunability of cold and ultracold molecular systems make them suitable to study a broad applications spectrum in diverse fields, from precision measurement and high-resolution spectroscopy \cite{Zelevinsky2008,DeMille2008} to ultracold chemistry \cite{Krems2005,Balakrishnan2016}, quantum information processing \cite{DeMille2002} and quantum computing~\cite{Micheli2006}.

Cold lattice gases, either atomic or molecular, are often proposed as quantum simulators mimicking, for instance, the Hubbard and Heisenberg \cite{Jaksch2005,Lewenstein2007,Rodriguez2010,Rodriguez2011} models, serving as a way to develop new experimental, theoretical and numerical tools. The production of ultracold ground-state molecules opens the path to derive effective many-body models for nonreactive molecules in tight traps like optical tweezers \cite{Kaufman2012,Thompson2013} and optical lattices \cite{Docaj2016}.  Arranging polar molecules in optical standing-wave trapping-potentials, new opportunities are expected due to the neighboring sites dipole-dipole interaction; in contrast, ground-state atomic lattice gases dominated by $1/R^6$ van der Waals forces exhibit only on-site interaction. Besides the long-range character of the dipolar interaction, it is interesting to explore its anisotropy which can be manipulated by applying external electric fields \cite{Rosario2007,Ni2010,Lameshko2013}.

The successful association of two pre-cooled atoms to form a molecule has been achieved using Feshbach resonances~\cite{Koehler2006} and photoassociation~\cite{Jones2006}. Magneto association performed by STIRAP has produced KRb molecules in their electronic and vibrational ground-state~\cite{Ni2008,Ospelkaus2008} and photoassociation close to a Feshbach resonance has produced LiCs molecules in their electronic and vibrational ground-state~\cite{Deiglmayr2008}.

It is of great interest to control the time-evolution among specific discrete quantum states, specifically, the array of long-range interacting polar molecules with rich internal structure form an excellent ground to comprehend complex systems relevant for quantum information processing and simulations~\cite{Covey2017}. For this porpoise, new theoretical and numerical tools are needed to understand the  many-body quantum systems such as rotational dynamics present in state of the art experiments.

Therefore, in the present work we analyze the coherent population transfer of $X^{1}\Sigma$ molecules due to the dipole-dipole interaction in their ground and first excited rotational levels as a function of an external electric field.  Pursuing this aim, the molecular gas is loaded onto an one-dimensional optical lattice and the gas is prepared in the strongly-correlated Mott-insulator regime.

Recently, the use of quantum information tools in several communities has been growth~\cite{Horodecki2009,Eisert2010}. Particularly, in the field of the many-body strongly-correlated quantum-systems~\cite{Amico2008,Islam2015}, these tools allow us to analyze nonlocal information and therefore a suitable mechanism to study quantum phases and phase transitions~\cite{mishmash2016,matsuura2016,He2017}. In a pure bipartite state, having measured only part of the system (a subsystem) it is possible to have information of the whole quantum state if there is entanglement among the parts~\cite{Calabrese2009}. In this paper, we use the von Neumann entropy associated to the reduced density matrix to provide the entanglement entropy of the measured subsystem.

This paper is organized as follows: in sect. \ref{model} we describe both the model of the molecular lattice system and the implemented numerical method. Section \ref{dynamics} is devoted to molecular excitations dynamics in presence of an external DC electric field. The entanglement generated throughout the evolution of the excitations is presented in sect. \ref{ent}. Finally, sect. \ref{summary} concludes.

\section{Model and methods} \label{model}

\subsection{Molecular Hamiltonian in the Mott-insulator regime} \label{hamiltonian}

We consider an optical lattice filled with $X^1 \Sigma$ LiCs mole\-cules which are in the electronic and vibrational ground states. The lattice holds one molecule per site at a distance of 400 nm, hence, tunneling between lattice sites is completely suppressed \cite{Herrera2010}. Considering this low lattice filling, the on-site Hubbard interaction can be neglected as well. Thus, the interaction between the molecules is solely determined by the long-range anisotropic dipole-dipole interaction, creating quasi-partic\-les responsible for excitation transfer along the lattice, so-called Frenkel excitons~\cite{Agranovich2008}. The only two states involved in the dynamics are the ground and the first excited rotational states. The Hamiltonian is given by~\cite{Herrera2010},

\begin{equation} \label{Excitation}
\hat H=\sum_{n=1}^{N_{\text{mol}}} \left(B_e\hat N_n^2-{\bold d}_n\cdot{\bold E}\right)+\frac{1}{2}\sum_{n=1}^{N_{\text{mol}}}\sum_{m\ne n}^{N_{\text{mol}}} \hat V_{dd}({\bold r}_n-{\bold r}_m),
\end{equation}

where ${\bold r}_n$ is the position of the $n$-th lattice site, $B_e$ is the rotational constant, ${\bold d}_n$ is the electric dipole operator of individual molecules, ${\bold E}$ is the constant electric field, $\hat{N}_n$ is the rotational angular momentum operator of the molecule at the site $n$ and $\hat{V}_{dd}$ is the dipole-dipole interaction between molecules in different lattice sites. $N_{mol}$ is the total number of molecules.

The local rotational Hamiltonian, $\sum_n B_e\hat N_n^2$, corresponds to the crystal energy operator (rotational energy) without intermolecular interactions. In absence of a static electric field, the rotational ground state of the given diatomic molecule has a rotational angular momentum $N=0$, and the first excited rotational state $N=1$ is triply degenerate. In the presence of a field  ${\bold E}$,  the Stark effect is accounted for, the degeneracy of $N=1$ is lifted and splits into two sub-energy levels corresponding to the projections, $M=0$ and $M=\pm 1$, into the electric field quantization axis.  For strong fields beyond a critical value (see reference \cite{Lincoln2009} for critical field theory), the system reduces to a two-level problem considering exclusively $M=0$ projections. In contrast, for a weak static electric field all four levels must be considered~\cite{Lincoln2009}, this is the scenario considered in the present work.

The electric field acts only on the internal coordinates of single molecules, and couples different local molecular states, this coupling depends on the electric field strength.

The dipole-dipole operator, $\hat V_{dd}({\bold r}_n-{\bold r}_m)$, involves two different sites and mixes two different molecular states. In fact this interaction determines the dynamics of the system and their features. In the Molecular Fixed Frame (MFF) this operator takes the form,

\begin{equation} \label{VddMFF}
\hat{V}_{dd}(\textbf{r})=\frac{1}{{r}^{3}}[\hat{d}_{n}\cdot\hat{d}_{m}-3(\hat{d}_{n}\cdot\hat{e})(\hat{e}\cdot\hat{d}_{m})].
\end{equation}

Here, $\textbf{r}={\bold r}_n-{\bold r}_m$ is the intermolecular vector, $\hat{e}=\frac{\bold{r}}{r}$ is the unit vector in the direction $\bold{r}$ and $\hat{d}_{m}, \hat{d}_{n}$ are the dipole moment operators of the molecules at the sites $m$ and $n$, respectively.

In order to study the interplay of the local interaction with $\textbf{E}$ and the non-local dipole-dipole interaction, it is needed to represent the operator (\ref{VddMFF}) in the Laboratory Fixed Frame (LFF),

\begin{equation} \label{VddLFF}
\hat{V}_{dd}(\textbf{r})=-\frac{\sqrt{6}}{{r}^{3}}\sum_{q=-2}^{2}{-1}^{q}C_{-q}^{2}(\textbf{r})[\hat{d}_{n}\otimes\hat{d}_{m}]_{q}^{(2)}.
\end{equation}

Here $C_{-q}^{2}(\textbf{r})=\sqrt{\frac{4\pi}{5}}Y^{(2)}_{-q}(\textbf{r})$ are the reduced spherical harmonics which describe the movement of the intermolecular axis in the LFF. The one-dimensional confinement of the lattice reduces this coefficient to a simple factor. The operators on the right hand side in Eq. (\ref{VddLFF}) contain the information about the strength of the involved molecular dipole moments. The term can be seen as a second rank tensor product of two tensors of rank one $\hat{T}_{i}^{1}(\hat{d}_{s})$ \cite{Zare} which represents the dipole moment operator of the molecule located in the site $s \in \{n,m\}$. As a result, we express $[\hat{d}_{n}\otimes\hat{d}_{m}]_{q}^{(2)}$ as $\sum_{p}\langle 1 p, 1 q-p | 2 q \rangle\hat{T}_{p}^{1}(\hat{d}_{n})\hat{T}_{q-p}^{1}(\hat{d}_{m})$.

Using the second quantization, the Hamiltonian (\ref{Excitation}) can be written following Ref.~\cite{Agranovich2008} as

\begin{equation} \label{finalH}
\hat{H}_{\text{exc}}=\sum_{n}E_{0}\hat{P}_{n,M}^{\dagger}\hat{P}_{n,M}+\sum_{n,m \neq n}\sum_{M,M'}F_{n,m}^{M,M'}\hat{P}_{n,M}^{\dagger}\hat{P}_{m,M'}.
\end{equation}

Here $E_{0}$ is the energy difference between the excited and ground rotational states, $\hat{P}_{n,M}^{\dagger}$ $(\hat{P}_{m,M'})$ creates (annihilates) a molecule in the site $n$ $(m)$ and level $M$ $(M')$. Specifically, the operators act as $\hat{P}_{n,M}^{\dagger} |g ,0\rangle _n = |e,M \rangle _n$ and $\hat{P}_{m,M'} |e,M' \rangle _m = |g \rangle _n$, where $|g,0 \rangle _n$ and $|e,M \rangle _n$ are the field-dressed ground and excited states, respectively. The coefficients $F_{n,m}^{M,M'}=\langle e_{n,M'} g_{m,0} | \hat{V}_{dd} | g_{n,0} e_{m,M} \rangle $ provide the molecular selection ru\-les, which determine the restricted phase space accessible at low-energies. Here, we use a DC electric field perpendicular to the molecular arrangement leading to the selection rule $\Delta (M_{n}+M_{m}) = 0, \pm 2$ involving molecules on sites $n$ and $m$~\cite{Herrera2010}. Finally, working with molecules separated distances greater than $100$ nm makes the intermolecular interaction considerably smaller compared to the rotational spectrum energy, so it is enough to account only for the nearest-neighbor interaction.

The Hamiltonian (\ref{finalH}) allows us to describe the behavior of the excitations as a function of the external electric field and study the coherent population transfer between the two involved rotational molecular levels.

We consider a lattice array with $L$ sites, such bipartite system can be separated by applying the Schmidt decomposition~\cite{Schmidt1907}. Each state can be written as $\ket{\Psi} = \sum_{\alpha_L =1}^{N_B} \lambda_\alpha \ket{\alpha_B} \otimes\ket{\alpha_{E}}$, where $\{\ket{\alpha_B}\}$ (on one subsystem with $N_B$ Hilbert dimension) and $\{\ket{\alpha_E}\}$ (on the other with $N_E$ Hilbert dimension) are two orthonormal basis sets belonging to the respective Hilbert subspace, each with at most $N_B$ basis elements, here $N_E>N_{B}$. The coefficients $\lambda_{\alpha}$ are non-negative real numbers satisfying $\sum_{\alpha}^{\chi}\lambda_{\alpha}^2=1$, known as the {\it Schmidt coefficients} and $\alpha$ runs in the reduced Hilbert subspace. In the numerical implementation the sum is truncated and only $\chi$ states are kept.

In a pure state there is only one density matrix eigenvalue. But the reduced density matrix for a subsystem $\hat\rho_B=\sum_\alpha\lambda_\alpha^2\ket{\alpha_B}\ket{\alpha_B}$ represents a mixed state where the entanglement entropy can be measured using the von Neumann entropy,

\begin{equation} \label{entropy}
S=-\sum_\alpha \lambda_\alpha^2 \ln (\lambda_\alpha^2).
\end{equation}

\subsection{Time evolution}

Numerical calculations reveal important aspects in the ultracold interacting molecular behavior.  In low dimensions, quantum fluctuations are enhanced and therefore many-body quantum systems, in general, are not susceptible to the standard mean field or perturbation expansion methods so one has to turn to non-perturbative techniques and numerical methods. Here, we use the Matrix Product State (MPS) ansatz~\cite{Vidal2003,Vidal2004,Verstraete2004,Verstraete-Porras2004} to efficiently approximate the quantum states. This formalism is highly recommended when the amount of entanglement in the system is limited, reducing the calculation complexity from the original exponential to an algebraic growth of the Hilbert space. Based on the MPS representation, we use the Time-Evolving Block Decimation (TEBD) algorithm~\cite{Vidal2004} to numerically analyze the molecular quantum-chain dynamics, paying special attention to population transfer and quantum entanglement. The TEBD takes advantage of the fact that the Hamiltonian can be written as the sum over {\it even} and {\it odd} sites and then the time evolution operator is approximated using the Trotter-Suzuki expansion formula~\cite{Suzuki1990}. Since the entanglement of the system grows with time, keeping the matrices of the same size means to loose information throughout the evolution. In order to avoid so, one has to adapt the value of the entanglement parameter ($\chi$) every time-step to keep the state as accurate as possible~\cite{Daley2004}.

The approximations discussed in sect.~\ref{hamiltonian} allow us to use this numerical technique for calculating the time evolution since the dipole-dipole interaction couples only nearest-neighbors~\cite{Rabl2007}. In the following, we consider chains with $L=16$ sites and initial entanglement parameter $\chi_{0}=8$ since the evolution starts from an array of disentangled individual molecules. Towards the end, the entanglement grows and we reach $\chi_{\text{final}}=100$ in average after 660 time-steps.

\section{Dynamics of the molecular excitations}\label{dynamics}

In this section, we first analyze the dynamics of the molecular excitations in the presence of a static and perpendicular electric field. Several scenarios are regarded depending on the field strength and the initial condition.

We work with LiCs alkali molecules, whose dipole moment is $d_0 = 5,523$ De~\cite{Aymar2005}, the rotational constant is $B_e = 5,816 $GHz and the critical electric field is $E_{\text{cri}} = 2100 $V/cm \cite{Lincoln2009}. Initially all lattice sites are prepared in the same state ($|\Psi_0\rangle$) which consists of either a single rotational level or a superposition of them. The choices are motivated by two ideas; first, the possibility to reach these single-molecule states after two-atoms association mechanisms. And second, the access to the strongly-correlated regime with dimers with tools developed in atomic physics field~\cite{Jaksch2005,Volz2006}.

As the system evolves the rotational excitations exchange and spread throughout the lattice following the selection rules. Hence, we locally measure the rotational state of each molecule, sum over all sites and divide by the total initial amount of molecules. Therefore, the figures presented in this section show the percent population of each single dressed state as a function of time. The time is given in unities of $(\mu^2/d^3)^{-1}$ which for our case is of the order of $\sim 9\times 10^4$ GHz$^{-1}$.

\begin{figure}
\psfrag{Population (\%)}{\hspace{-0.3cm}\footnotesize Population $(\%)$}
\psfrag{Time}{\hspace{-0.6cm}\footnotesize Time $(\mu^2/d^3)^{-1}$}
\psfrag{0,0}{{\hspace{-0.3cm}\scriptsize  $|0,0\rangle$}}
\psfrag{1,-1}{{\hspace{-0.3cm}\scriptsize  $|1,\text{-}1\rangle$}}
\psfrag{1,0}{{\hspace{-0.3cm}\scriptsize $|1,0\rangle$}}
\psfrag{1,1}{{\hspace{-0.3cm}\scriptsize $|1,1\rangle$}}
\psfrag{(a)}{{\hspace{-1.1cm}\scriptsize (a) $|0,0\rangle$}}
\psfrag{(b)}{{\hspace{-1.1cm}\scriptsize (b) $|1,0\rangle$}}
\psfrag{(c)}{{\hspace{-1.1cm}\scriptsize (c) $|1,1\rangle$}}
\psfrag{(d)}{{\hspace{-1.1cm}\scriptsize (d) $|\Psi_1\rangle$}}
\psfrag{(e)}{{\hspace{-1.1cm}\scriptsize (e) $|\Psi_2\rangle$}}
\psfrag{(f)}{{\hspace{-1.1cm}\scriptsize (f) $|\Psi_3\rangle$}}
\includegraphics[width=1\linewidth]{mult_5000.eps}
\caption{Coherent population transfer among the rotational states in presence of a weak DC electric field of $E=50$ V/cm. All molecules are prepared in the same Fock state (a) $|0,0\rangle$, (b) $|1,0\rangle$ and (c) $|1,1\rangle$ and the superposition (d) $|\Psi_1\rangle=\frac{1}{\sqrt{2}}\big(|1,\text{-}1 \rangle + |1,1 \rangle\big)$, (e) $|\Psi_2\rangle=\tfrac{1}{\sqrt{19}}\big(3|1,\text{-}1\rangle + |1,0 \rangle+ 3|1,1 \rangle\big)$ and (f) $|\Psi_3\rangle=\frac{1}{\sqrt{3}}\big(|1,\text{-}1 \rangle + | 1,0 \rangle + | 1,1 \rangle\big)$. See color figures online.\label{E50}}
\end{figure}

\subsection {Weak field: Effective two- and three-level dynamics}

Let us start by analyzing the case where a weak field is applied, namely $E=50 $ V/cm.
For a very low external field, the coupling between $| 0,0 \rangle$ and $| 1,0 \rangle$ is weak. This physical fact is presented in Fig.~\ref{E50}, starting from the $| 0,0 \rangle$ state (see Fig. \ref{E50}(a)), the population rapidly starts to decrease and head towards the two degenerate states $|1,-1 \rangle$ and $| 1,1 \rangle$, while $|1,0 \rangle$ state remains unchanged up to $t=0.5$ and slightly, almost imperceptible, rises its population.
On the contrary, when all the population starts in the $|1,0 \rangle$ state, the excluded states from the exchange are $| 1,-1 \rangle$ and $| 1,1 \rangle$ (see Fig. \ref{E50}(b)). The system behaves as a {\it two-level system}, besides the initial population does not decrease dramatically as in  \ref{E50}(a), but slowly go to the coupled state $|0,0 \rangle$. Initiating from $|1,1 \rangle$  the population first go to the ground state due to the selection rules and later on the state $| 1,-1 \rangle$ begins to be populated. It can be seen how $| 1,-1 \rangle$ takes a while before enters into the dynamics and afterwards $| 0,0 \rangle$, $| 1,-1 \rangle$ and $| 1,1 \rangle$ coherently oscillate.

Otherwise, Figs.~\ref{E50}(d), (e) and (f) show different initial population configurations, the evolutions start from linear superpositions of $| 1,-1 \rangle$, $| 1,0 \rangle$ and $| 1,1 \rangle$ states with different contributions.  For these three cases the systems behave effectively as {\it three-level systems} involving only the states: $\{|0,0\rangle, |1,-1\rangle, |1,1\rangle\}$, excluding the contribution of the $|1,0\rangle$ state, which keeps the population practically constant for the whole evolution. The same behavior applies for Figs.~\ref{E50}(a) and (c).

Hence, this low DC field strength allows us, depending on the initial condition, to switch from a four-level system to a two- or a three-level one.

\subsection {Moderate field}
\begin{figure}[t
\psfrag{Population (\%)}{\hspace{-0.4cm} \footnotesize Population $(\%)$}
\psfrag{Time}{\hspace{-0.6cm}\footnotesize Time $(\mu^2/d^3)^{-1}$}
\psfrag{0,0}{{\hspace{-0.3cm}\scriptsize  $|0,0\rangle$}}
\psfrag{1,-1}{{\hspace{-0.3cm}\scriptsize  $|1,\text{-}1\rangle$}}
\psfrag{1,0}{{\hspace{-0.3cm}\scriptsize $|1,0\rangle$}}
\psfrag{1,1}{{\hspace{-0.3cm}\scriptsize $|1,1\rangle$}}
\psfrag{(a)}{{\hspace{-1.1cm}\scriptsize (a) $|0,0\rangle$}}
\psfrag{(b)}{{\hspace{-1.1cm}\scriptsize (b) $|1,0\rangle$}}
\psfrag{(c)}{{\hspace{-1.1cm}\scriptsize (c) $|1,1\rangle$}}
\psfrag{(d)}{{\hspace{-1.1cm}\scriptsize (d) $|\Psi_1\rangle$}}
\psfrag{(e)}{{\hspace{-1.1cm}\scriptsize (e) $|\Psi_2\rangle$}}
\psfrag{(f)}{{\hspace{-1.1cm}\scriptsize (f) $|\Psi_3\rangle$}}
\includegraphics[width=1\linewidth]{mult_20000.eps}
\caption{Coherent evolution starting all molecules from the initial conditions defined in Fig.~\ref{E50} and using a moderate DC electric field strength of $E=200$ V/cm. See color figures online.\label{E200}}
\end{figure}

As the field strength rises, all four states are involved in the dynamics. Here we use $E=200$ V/cm and the coherent collective evolutions are shown in Fig. \ref{E200}. Starting from the state $|0,0 \rangle$ (Fig.~\ref{E200}(a)) or state $|1,0 \rangle$ (Fig.~\ref{E200}(b)), the difference directly regards the selection rules. In the case of \ref{E200}(a) the excitations migrate directly to $N=1$, whereas in the case of \ref{E200}(b) the excitations must first go to the ground state and thereafter  populate the $M=\pm 1$ states. Hence, the states $|1, \pm 1 \rangle$ remain unpopulated for a longer period of time.

The coupling between $M=0$ projections gets stronger than the latter case ($E=50$ V/cm), so in Fig.~\ref{E200}(a) the population of $|1,0 \rangle$ state heightens up to $20\%$ by the end of the time evolution, and the smooth coherent oscillations between $| 0,0 \rangle$ and $| 1,\pm 1 \rangle$ observed in Fig.~\ref{E50}(a) are weakened even more for $E=200$ V/cm. On top of that, after $t=1.5$ the populations in $|1, \pm 1 \rangle$  level off to a stationary regime getting 25\% of the total population each. From this point on the dynamics reach again an effective {\it two-level like system}.
In Fig.~\ref{E200}(b) the interaction between ground $N=0$ and exited states $N=1$ increases and the exchanges is higher than the case with lower electric field \ref{E50}(b), and slight oscillations among the rotational levels appear.

When the $|1,1 \rangle$ state is taken to be the initial state (similar analysis applies for $|1,-1 \rangle$), see Fig. \ref{E200}(c), the population rapidly drops into $|0,0 \rangle$ state, then goes into $|1,-1 \rangle$ and finally the $|1,0 \rangle$ enters in the dynamics. The field dressed states with $M = \pm 1$ have the same energy, hence, they have the same probability to be populated, being the reason why $|1,-1 \rangle$ state comes first in the dynamics than $|1,0 \rangle$ state. High-amplitude coherent oscillations are presented for this initial condition and the considered field strengths.

For the case shown in Fig.~\ref{E200}(d), the evolution starts from a superposition state with equal populations (symmetric Bell-like state) using $|1,\pm 1 \rangle$ instead of a Fock state. As it is observed, the rotational ground state becomes quickly populated while $|1,0 \rangle$ rises late and slowly in the dynamics. The latter levels out reaching a steady regime whilst the populations of the other three states oscillate coherently, regardless of the initial Fock state (c) or Bell superposition (d).

As long as $|1,\pm 1 \rangle$ dressed states are initially populated, the $M=0$ projection of the $N=1$ state behaves in the same manner matching actually their time-scales. This observation still applies even when the $|1,0 \rangle$ state is allowed to be populated with almost 10$\%$ of the total population, see Fig. \ref{E200}(e) (calculations where performed also with 20$\%$ of the population delivering the same result) and the behavior for this state persists: remaining constant at the initial population for a while, and then rising slowly until it levels out.

Finally we explore the case when all states in the rotational excited state $N=1$ have the same initial weight achieving the evolution presented in Fig. \ref{E200}(f). In this case, the rotational level $|1,0\rangle$ has the largest probability of being populated as the evolution takes place, while the populations in the ground state and projections $M=\pm 1$ tend to a quasi-steady regime sharing the same probability percentage.
Between $t=0.5$ and $t=1$, the ground state reaches its maximum, almost when the population in $| 1,0 \rangle$ stars to rise, see Figs.~\ref{E200}(c), (d) and (f). Therefore, for the ground state to get populated a elapsed time around $(\mu^2/d^3)^{-1}$ is needed, afterwards a redistribution to the other leves is performed according to selection rules and weight probabilities.

\subsection {Strong field}
\begin{figure}[t]
\psfrag{Population (\%)}{\hspace{-0.3cm}\footnotesize Population $(\%)$}
\psfrag{Time}{\hspace{-0.6cm}\footnotesize Time $(\mu^2/d^3)^{-1}$}
\psfrag{0,0}{{\hspace{-0.3cm}\scriptsize  $|0,0\rangle$}}
\psfrag{1,-1}{{\hspace{-0.3cm}\scriptsize  $|1,\text{-}1\rangle$}}
\psfrag{1,0}{{\hspace{-0.3cm}\scriptsize $|1,0\rangle$}}
\psfrag{1,1}{{\hspace{-0.3cm}\scriptsize $|1,1\rangle$}}
\psfrag{(a)}{{\hspace{-1.1cm}\scriptsize (a) $|0,0\rangle$}}
\psfrag{(b)}{{\hspace{-1.1cm}\scriptsize (b) $|1,0\rangle$}}
\psfrag{(c)}{{\hspace{-1.1cm}\scriptsize (c) $|1,1\rangle$}}
\psfrag{(d)}{{\hspace{-1.1cm}\scriptsize (d) $|\Psi_1\rangle$}}
\psfrag{(e)}{{\hspace{-1.1cm}\scriptsize (e) $|\Psi_2\rangle$}}
\psfrag{(f)}{{\hspace{-1.1cm}\scriptsize (f) $|\Psi_3\rangle$}}
\includegraphics[width=1\linewidth]{mult_100000.eps}
\caption{Rotational levels dynamics with all molecules prepared in the initial conditions defined in Fig.~\ref{E50} using a strong DC electric field of $E=1000$ V/cm. See color figures online.\label{E1000}}
\end{figure}

Let us now increase the electric field to $E=1000$ V/cm which is still small compared to the rotational molecular energy. In this scenario, the coupling between the states with $M=0$ projection gets stronger, therefore the dynamics tend to be dominated by these two states, while the $M=\pm 1$ projections behave smoothly and the evolutions for several initial states are presented in Fig.~\ref{E1000}.

Let us consider two initial cases, when all molecules are prepared in the rotational ground state (Fig. \ref{E1000}(a)) and in the first rotational excited state with $M=0$ (Fig. \ref{E1000}(b)). In both cases the initial state drops dramatically its population  as much as the coupled state with the same projection rises and the populations oscillate strongly between the $m=0$ projections. Thereafter, these two states exchange coherently most of the chain population while the  $|1, \pm 1 \rangle$ dressed states starting from zero have an upward trend, with a slight decrease, stabilizing towards the end. For the evolution \ref{E1000}(b) $M=\pm 1$ states remain unpopulated for a longer period of time climbing up when the populations in $|0,0\rangle$ and $|1,0\rangle$ match each other, leveling out although swinging slightly until the end of the evolution.

Taking $|1,1 \rangle$ as the initial state leads to a different time evolution (see Fig.~\ref{E1000}(c)). According to the selection rules, molecules go first to the ground state and since the $M=0$-coupling becomes stronger due to the hight field, the $|1,0 \rangle$ state has a larger probability to be populated. Hence, in contrast of the previous case (Fig.~\ref{E200}(c)), the level population order changes and the ground state enters first in the dynamic, then the $| 1,0 \rangle$ state takes advantage over $| 1,-1 \rangle$ that comes latter. Once all levels get populated the states with $m=0$ oscillate coherently, while the $|1, \pm 1 \rangle$ states are struggling to equate their populations. The maximum of  $| 1,- 1 \rangle$ is reached approximately at $t\sim 2$.

When launching the evolution with each molecule prepared in a Bell state made up of $|1,\pm 1\rangle$ (see Fig.~\ref{E1000}(d)), the population in both levels decreases continuously during the first time unity and remains stable around 30\% throughout the rest of the evolution. Meanwhile, the $|0,0 \rangle$ and $|1,0 \rangle$ states behave in the same manner as in Fig.~\ref{E1000}(c). Modifying the initial state by increasing the population in $|1,0 \rangle$ does not change the general demeanor of the evolution of the $M=\pm 1$ projections as shown in Fig.~\ref{E1000}(e). This is true even if the initial population in $|1,0 \rangle$ equals $M=\pm 1$ as observed in Fig.~\ref{E1000}(f). In the latter, we analyze the initial state with the molecules in a linear combination of equally probable $N=1$ states. Here, one can see that after decreasing the population in $M=\pm 1$ level almost to a constant value whilst both $M=0$ projections oscillate coherently damping down towards the end of the evolution.

Moderate and strong field strengths lead, regardless the initial conditions, the evolutions to a population distribution between $20 \%$ and $40 \%$. Hence, for longer period of times, out of our calculations scope, we infer that the population  distribution in the steady state is expected to land in the same range.

\section{Entanglement evolution.} \label{ent}

\begin{figure}[t]
\psfrag{Time}{\hspace{-0.6cm}\footnotesize Time $(\mu^2/d^3)^{-1}$}
\psfrag{0s}{{\hspace{-0.21cm}\scriptsize $|0,0\rangle$}}
\psfrag{0t}{{\hspace{-0.228cm}\scriptsize $|1,0\rangle$}}
\psfrag{1}{{\hspace{-0.31cm}\scriptsize $|1,1\rangle$}}
\psfrag{111}{{\hspace{-0.1cm}\scriptsize $|\Psi_1\rangle$}}
\psfrag{313}{{\hspace{-0.1cm}\scriptsize $|\Psi_2\rangle$}}
\psfrag{Entanglement}{\hspace{-0.6cm}\scriptsize Entanglement entropy}
\includegraphics[width=\linewidth]{ent_ver.eps}
\caption{Von Neumann entanglement entropy generated in the dynamics presented in Figs.~\ref{E50}, \ref{E200} and \ref{E1000} using as initial states: the Fock states $\{|0,0\rangle,|1,0\rangle,|1,1\rangle\}$ and the superpositions $\{|\Psi_1\rangle,|\Psi_2\rangle\}$ with $|\Psi_1\rangle=\tfrac{1}{\sqrt{3}}\big(|1,\text{-}1 \rangle + | 1,0 \rangle+ | 1,1 \rangle\big)$ and $|\Psi_2\rangle=\tfrac{1}{\sqrt{19}}\big(3| 1,\text{-}1 \rangle + |1,0 \rangle + 3|1,1 \rangle\big)$. Color figures online.\label{Ent}}
\end{figure}

In Sec.~\ref{dynamics} the evolution of the rotational levels has been presented depending on the initial state. Now, we focus on how the von Neumann entanglement entropy of the system grows as a function of time for three specific DC field strengths and several initial conditions. The corresponding results are shown in Fig.~\ref{Ent}.

At $E=50$ V/cm, see Fig.~\ref{Ent}(a), the system develops the highest but also the lowest entanglement observed in our numerical calculations.
 Starting all molecules in the rotational ground-state turns out to be the initial condition that most entangles the system (lila curve), this is due to the fact that this state interacts with all the excited states, backed up by the selection rules established by the dipole-dipole interaction. Let us recall that for this electric field strength the state $|1,0\rangle$ is excluded while the populations in $M=\pm 1$ states increase until the three involved states level out reaching a quasi-stationary regime. Therefore, it is expected that initiating from $|1,0\rangle$ (green curve) the evolution leads to the most disentangled dynamic and this is exactly what we observed.
 When the majority contribution to the initial state comes from $ M= \pm 1$ states (red and blue curves), independently if it is a Fock or superposition state, due to the selection rules these levels have access only to the ground state, thus the system develops similar entanglement throughout the evolution.
 Even when the initial state, $|\Psi_2\rangle=\tfrac{1}{\sqrt{19}}\big(3| 1,\text{-}1 \rangle + |1,0 \rangle + 3|1,1 \rangle)$, which has a partial contribution of $| 1,0 \rangle$, the greatest contribution  to the dynamics comes from $M = \pm 1$ states, which determine the entanglement evolution.
A different scenario is presented when launching the evolution from an equally probable superposition of the excited rotational states (orange curve), where the entanglement grows steadily changing the slope to a lower rate at two time-units. This is because such superposition connects exclusively to the ground state by the selection rules.

As the field strength is increased to $E=200$ V/cm, see Fig.~\ref{Ent}(b), all evolutions rapidly entangle except the one when taking $|1,0\rangle$ as initial condition (green curve). In the dynamic this initial state mostly interacts with $ | 0,0 \rangle$ state, while the interaction with $m=\pm1$ is through of  $| 0,0 \rangle$. Notice that initiating from $| 1,1 \rangle$ and the superposition $|\Psi_1\rangle=\tfrac{1}{\sqrt{3}}\big(|1,\text{-}1 \rangle + | 1,0 \rangle + | 1,1 \rangle\big)$  lead to two different dynamics (Figs.~\ref{E200}(c) and (f)), but similar entanglement evolution, since both states transfer their population directly to the ground state. Let us point out that the coupling between $| 1,0 \rangle$ and $| 0,0 \rangle$ gets enhanced as the external electric field increases, hence it differs from the interaction among the levels $| 1, \pm 1 \rangle$ and $| 0, 0 \rangle$.

 Finally in the case of $E=1000$ V/cm, see Fig.~\ref{Ent}(c), most of the curves coincide except for the evolution of the superposition with equally probable $N=1$ states whose evolution does not match the rest although keeps the same behavior below the others. For this electric field the coupling between $M=0$ states gets even stronger, and $M=\pm 1$ states gets weaker, then for superposition states (blue and orange curves) that include equally probable $M=\pm1$ states, their entanglement evolution slow down.

 Overall, it is observed that the weaker the field the higher the entanglement rises for the ground state, this is due to its connection with all $N=1$ states. As the electric field increases, the coupling between $M=0$ levels is favored over $M=\pm1$ states and the amount of states involved throughout the evolution decreases and the entanglement in the chain decreases as well.

Let us note that at some point all curves must stop the growing tendency and level out which is not presented in our calculations. This is due to the fact that the entanglement for finite size systems do not grow indefinitely rather reach a logarithmic behavior in 1D, and our calculations show only an early evolution of the molecular chain.

\section{Conclusions} \label{summary}

We studied a LiCs molecular gas loaded on an 1D optical lattice prepared in the strongly-correlated Mott-insulator regime with one molecule per lattice site.  The molecules are in the lowest electronic and vibrational state and the coherent population transfer between the internal rotational states as a function of an external DC electric field is analyzed.  Due to the large intermolecular lattice distance and low filling, the hopping between nearest-neighbors and the on-site Hubbard interaction are suppressed, hence, the dipole-dipole interaction in the near\-est-neighbor approximation solely governs the dynamics of the excitations. For field strengths lower than the critical field the full set of rotational levels must be taken into account and the numerical simulations, here reported, confirmed this fact. Coherent transfer of population among the field-dressed states along the lattice are shown for several electric field intensities and initial conditions. Our numerical simulations show that weak fields have the potential to peak either a two- or three-level system out of the original four-level one, therefore, the field acts as field-selector for the molecular system, simplifying the complexity for both numerical and analytical calculations. Even more, the excluded states can be thought of as dark states and can be used to store unperturbed information throughout the dynamics regarding the initial condition.

Numerical simulations, performed with three field stre\-ngths, showed that molecular excitations become entangled as the dynamics occur for most initial conditions. The higher the field intensity, the more certainly the system will become entangled, regardless of the initial condition. Moreover, a sharp and monotonic growth of the entanglement was in general obtained for early time evolutions, reason why these type of novel molecular systems are suitable to be proposed for new quantum protocols.

\section*{Acknowledgements}

K.R. thanks Ian Picken and Christiane P. Koch for carefully reading the manuscript and their valuable comments. This work has been supported by Universidad del Valle under the internal project 7967. The authors acknowledge the support from the Colombian Science, Technology and Innovation Fundation--COLCIENCIAS ``Francisco Jos{\'e} de Caldas'' under project 110\-665842793 (contract FP-005-2015). We also thank the Center of Excellence for Novel Materials (CENM) at Universidad del Valle for financial support for the research group.

\bibliographystyle{ieeetr}

redpajama_set_name: RedPajamaArXiv


text:
South Carolina Highway 341 (SC 341) is a  state highway that travels through the northeastern part of the U.S. state of South Carolina. It travels diagonally from the southeast to northwest (but signed as north–south) from SC 41 and SC 51 in Johnsonville to U.S. Route 521 Business (US 521 Bus.) and US 601 Bus. in Kershaw.

Route description

SC 341 begins at the intersection of Georgetown Highway (SC 41/SC 51) and Broadway Street in central Johnsonville. SC 341 heads west along East Broadway Street through the town's business district. At a railroad crossing, the road name changes to West Broadway and continues past additional businesses. After exiting the city limits, the road's name changes to Lake City Highway but passes some other small businesses, houses, and churches. Further from the city, more farmland begins to surround the roadway. At the community of Prospect Corner, the road name becomes Johnsonville Highway but the highway's surroundings is generally the same. Before entering Lake City, SC 341 passes in front of some commercial and industrial warehouses. Upon entering Lake City, SC 341 travels through it on Main Street, first through a residential neighborhood with an elementary school located along it. At Church Street, US 378 Bus. joins SC 341 on a concurrency through the city. Commercial businesses line both sides of the street upon the beginning of the concurrency. In the center of the city, the two routes intersect US 52. After passing Matthews Road west of the downtown area, the road enters another residential neighborhood but then heads into an area with more farmlands. Upon its intersection with US 378, US 378 Bus. ends and SC 341 has a short  concurrency with US 378.

Continuing northwest from US 378, SC 341 travels along Olanta Highway. At the community of Camp Branch, the highway intersects SC 403 at its southern terminus. SC 341 continues through mostly farmland and afterwards intersects SC 541 at its southern terminus. Shortly thereafter, the road enters the town of Olanta on Hampton Street. Within the town, the road passes in front of many houses, some small businesses, a library, and a school. In the center of town, SC 341 intersects US 301 then begins to curve to the north. The road heads through a mix of farmland and woods north of Olanta where it exits Florence County and enters Sumter County. After an intersection with SC 53, SC 341 begins to curve to the northwest and has an interchange with Interstate 95 (I-95) at its exit 146. After entering Lee County, SC 341 reaches the community of South Lynchburg where it intersects SC 58 at its western terminus. Soon after, it enters the town of Lynchburg where in the center, it has an intersection with US 76. It continues northwest along Lynchburg Highway through farmland and intersects US 401 about  northeast of Elliott. At the settlement of Wisacky, SC 341 reaches the northern terminus of SC 527.

Now approaching the city of Bishopville, the highway passes in front of Lee Central High School and has an interchange with I-20 at exit 120. North of the Interstate, it passes the entrance to Lee Correctional Institution. Upon entering the city of Bishopville, it takes the name of East Church Street where it heads through a residential neighborhood then transitions to the business district. In the town center, SC 341 reaches Main Street which carries US 15 through the city. SC 341 and SC 34 (which came into this intersection from West Church Street) form a concurrency with US 15 and head northeast through Bishopville. Outside the city limits, SC 341 breaks off the US 15/SC 34 concurrency and heads north-northwest through farmland and woods on Bethune Highway. The highway enters Kershaw County and soon afterwards heads into the town of Bethune. It travels through the town on Main Street where it intersects US 1 and SC 346. Outside of Bethune, it heads through the pine forests of the Sandhills region. After crossing the Little Lynches River, the road enters Lancaster County. It then enters the town of Kershaw on South Minor Street. At Marion Street, SC 341 turns to the west while SC 157 begins and travels east on Marion Street and SC 157 Truck / SC 341 Truck head north on North Minor Street. The highway heads into the center of town where SC 341 ends at Hampton Street which carries US 521 Bus./US 601 Bus. Marion Street continues west to intersect US 501 and US 601 at Matson Street.

Major intersections

Kershaw truck route

South Carolina Highway 341 Truck (SC 341 Truck) is a  truck route for trucks traveling through Kershaw. For its entire length, it is concurrent with SC 157 Truck. Its southern terminus is at the intersection of Marion Street and Minor Street which is also where SR 341 has its southern terminus and SR 157 travels through. As trucks are banned on East Marion Street through Kershaw, the truck route bypasses the downtown area. The highways head north along North Minor Street for three blocks to Hilton Street, which carries U.S. Route 601 (US 601). The northbound directions of the two truck routes travel along southbound US 601 forming a wrong-way concurrency. At North Hampton Street, the highways reach U.S. Route 521 Business (US 521 Bus.) and US 601 Bus. Continuing west along Hilton Street, they cross a railroad track before the truck routes terminate at US 521 (North Matson Street). US 601 from this point continues south along US 521 around the west side of the town.

See also

References

External links

SC 341 at Virginia Highways' South Carolina Highways Annex
Former SC 341 ALT at Virginia Highways' South Carolina Highways Annex

341
Transportation in Florence County, South Carolina
Transportation in Sumter County, South Carolina
Transportation in Lee County, South Carolina
Transportation in Kershaw County, South Carolina
Transportation in Lancaster County, South Carolina

redpajama_set_name: RedPajamaWikipedia