
Citation: Zareen Gul, Muhammad Younas Khan Barozai, Muhammad Din. In-silico based identification and functional analyses of miRNAs and their targets in Cowpea (Vigna unguiculata L.)[J]. AIMS Genetics, 2017, 4(2): 138-165. doi: 10.3934/genet.2017.2.138
[1] | José M. Campos-Salazar, Roya Rafiezadeh, Felipe Santander, Juan L. Aguayo-Lazcano, Nicolás Kunakov . Comprehensive GSSA and D-Q frame dynamic modeling of dual-active-bridge DC-DC converters. AIMS Electronics and Electrical Engineering, 2025, 9(3): 288-313. doi: 10.3934/electreng.2025014 |
[2] | Syeda Nadiah Fatima Nahri, Shengzhi Du, Barend J. van Wyk, Oluwaseun Kayode Ajayi . A comparative study on time-delay estimation for time-delay nonlinear system control. AIMS Electronics and Electrical Engineering, 2025, 9(3): 314-338. doi: 10.3934/electreng.2025015 |
[3] | Vengadesan Alagapuri, Ashok Bakkiyaraj Radhakrishnan, S. Sakthivel Padaiyatchi . Optimization of position and rating of shunt and series connected FACTS devices for transmission congestion management in deregulated power networks. AIMS Electronics and Electrical Engineering, 2024, 8(2): 165-186. doi: 10.3934/electreng.2024007 |
[4] | Muamer M. Shebani, M. Tariq Iqbal, John E. Quaicoe . Comparison between alternative droop control strategy, modified droop method and control algorithm technique for parallel-connected converters. AIMS Electronics and Electrical Engineering, 2021, 5(1): 1-23. doi: 10.3934/electreng.2021001 |
[5] | Yanming Wu, Zelun Wang, Guanglei Meng, Jinguo Liu . Neural networks-based event-triggered consensus control for nonlinear multiagent systems with communication link faults and DoS attacks. AIMS Electronics and Electrical Engineering, 2024, 8(3): 332-349. doi: 10.3934/electreng.2024015 |
[6] | Mohammed Mecifi, Abdelmadjid Boumediene, Djamila Boubekeur . Fuzzy sliding mode control for trajectory tracking of an electric powered wheelchair. AIMS Electronics and Electrical Engineering, 2021, 5(2): 176-193. doi: 10.3934/electreng.2021010 |
[7] | Mohamed M. Albarghot, Mohamed T. Iqbal, Kevin Pope, Luc Rolland . Dynamic modeling and simulation of the MUN Explorer autonomous underwater vehicle with a fuel cell system. AIMS Electronics and Electrical Engineering, 2020, 4(1): 114-131. doi: 10.3934/ElectrEng.2020.1.114 |
[8] | Tadele A. Abose, Thomas O. Olwal, Muna M. Mohammed, Murad R. Hassen . Performance analysis of insertion loss incorporated hybrid precoding for massive MIMO. AIMS Electronics and Electrical Engineering, 2024, 8(2): 187-210. doi: 10.3934/electreng.2024008 |
[9] | José M. Campos-Salazar, Juan L. Aguayo-Lazcano, Roya Rafiezadeh . Non-ideal two-level battery charger—modeling and simulation. AIMS Electronics and Electrical Engineering, 2025, 9(1): 60-80. doi: 10.3934/electreng.2025004 |
[10] | Lawrence O. Aghenta, M. Tariq Iqbal . Design and implementation of a low-cost, open source IoT-based SCADA system using ESP32 with OLED, ThingsBoard and MQTT protocol. AIMS Electronics and Electrical Engineering, 2020, 4(1): 57-86. doi: 10.3934/ElectrEng.2020.1.57 |
Radiocarbon (14C) dating represents a key chronological method that has been widely used for determining the age of geological and archaeological samples of the last ca. 55,000 years [1]. The physical basis of this method is that organisms assimilate 14C, directly or indirectly, from the atmosphere during their life span. However, their death cuts off the linkage to the contemporaneous global carbon cycle, which leaves the remainder 14C atoms inside their bodies to decay without replenishment with a half-life of ca. 5730 years. Therefore, by the law of radioactive decay, the radiocarbon age of a sample can be determined uniquely by measuring the residual 14C radioactivity [2], providing that the atmospheric 14C level is assumed to remain constant over time. However, owing to the complex global carbon cycle and the changes in the intensity of the solar and earth's magnetic field [3], the atmospheric 14C level fluctuated remarkably [4], exerting a profound impact on the radiocarbon dating method [5]. Therefore, the laboratory radiocarbon ages must be calibrated using a calibration curve [6], which is a key step for converting radiocarbon ages to calendar ages'.
Due to the non-monotonic nature of the calibration curves, the map from the radiocarbon year to the calendar year is not one to one, thereby leading to non-unique solutions usually expressed in terms of a probability density function with one or several highest posterior density (HPD) regions [7,8]. It appears that uncertainty is a nature of the calibrated age, which cannot be eliminated regardless of how small the laboratory error is [5]. Sometimes the uncertainty of the calibrated age could be extremely large if the radiocarbon age hits a plateau in the calibration curve [9,10]. Many attempts have been made to reduce the uncertainty of radiocarbon age calibration. Bayesian methods represent a promising approach to the solution of this problem [11,12,13], which takes advantage of the stratigraphical context or temporal order of the radiocarbon ages as a constraint [14]. Specifically, if the radiocarbon ages are assumed to come from a depositional sequence, they must follow a chronological order [15,16]. By applying the constraint of temporal order on the calibration, the uncertainty of the calibrated age can be greatly narrowed down [8].
Bayesian methods have been increasingly used for modeling a suite of 14C ages subject to temporal order constraints since 1990s. Several Bayesian frameworks have been developed and in use. Of particular note are the standalone package DateLab [17] and the online OxCal [11] and BCal [18] tools. DateLab provides a simple graphical user interface, while OxCal and BCal allow registered users to access through a worldwide web (www) browser for data entry and workflow setup. However, further analyses of the modeled ages require the random samples from the original Markov Chain Monte Carlo (MCMC) output. For example, if two ages need to be compared, the best way is to calculate the difference between their posterior ages and expresses it in terms of a probability density function with the HPD regions. Also, inexperienced users may need not only a more flexible and user-friendly environment for data entry and definition of their project-specific problems, but also a powerful post-processing tool for analyzing and visualizing the results. Towards this end, a hierarchical Bayesian framework for 14C age modeling subject to temporal order constraints is described and a Matlab package is presented here. The package is an extension to the previous version [19] by adding several new features such as the overlapping boundary relationship that deal with the coexistence and subsequent replacement of one period by another.
Considering a chronological sequence consisting of M periods, where "period" is defined as a stage of deposition or human occupation characterized by a distinct geological or archaeological feature. Similar to the "phase" in OxCal, each period contains a group of at least one radiocarbon age. For the ith period of the chronological sequence, where i=1,2,⋯,M denoting the number of the lowest (oldest) to the uppermost (most recent) period, Let ri=[ri,1,ri,2,⋯,ri,Ni], σσi=[σi,1,σi,2,⋯,σi,Ni], θθi=[θi,1,θi,2,⋯,θi,Ni], αi, and βi denote the Ni 14C ages, the associated one standard deviation of the Ni 14C ages, the corresponding calendar ages of the Ni 14C ages, and the calendar age of the early and late boundary of the ith period, respectively (Figure 1). Note that the division of periods in the sequence is subjective, but mainly depends on the sedimentary facies in geology or tool assemblage in archaeology or the scientific questions to be answered. Here a hierarchical Bayesian framework of 14C age modeling with a minimum level of structural complexity is described (Figure 2), within which θi, αi, and βi, where i=1,2,⋯,M, are treated as model parameters, which will be inferred using the laboratory radiocarbon ages. Due to lack of knowledge about these parameters, vague (noninformative) priors are preferred to use. For the sake of mathematical expression, the "BP" scale is used for model description, while the algorithms were implemented in both "BP" and "BC/AD" scales. The conversion from the "BP" to the "BCE|CE" scale is BCE|CE = 1950 – BP.
For the Ni radiocarbon ages in the ith period, where i=1,2,⋯,M, the corresponding calendar ages, say θi, are assumed to follow the uniform distribution supported on the interval of [αi,βi].
f(θθi|αi,βi)=1(αi−βi)NiI(θθi|αi,βi), | (1) |
where I(θθi|αi,βi) is an indicator function describing the temporal order of the calendar ages in the ith period. Here, three types of temporal order are considered: (1) "unordered", which means that the ages are not necessarily to be chronologically ordered in a period. An example is that ages were obtained from different archaeological sites of the same cultural period. As such, the indicator function can be defined as:
I(θθi|αi,βi)={1,αi>all(θθi)>βi0,otherwise, | (2) |
(2) "coeval", which means that the ages are assumed to be the same in a period. An example is that ages were obtained from the same depth of a stratigraphical unit or the same archaeological feature of a cultural period and thus their calendar ages must be the same. In this case, the indicator function can be defined as:
I(θθi|αi,βi)={1,αi>θi,1=θi,2=⋯=θi,Ni>βi0,otherwise, | (3) |
and (3) "ordered", which means that the ages follow a chronological order. An example is that ages were obtained from different depths of a stratigraphical unit. Therefore, the indicator function can be expressed as:
I(θθi|αi,βi)={1,αi>θi,1>θi,2>⋯>θi,Ni>βi0,otherwise, | (4) |
where " > " denotes "older than" and "=" denotes "same as".
The calendar age of the early and late boundaries of the ith period, say αi and βi, where i=1,2,⋯,M, is assumed to follow the uniform distribution supported on the interval [A, B], where A and B are two floating parameters defining the beginning and the end of the chronological sequence in calendar years BP, respectively. Therefore, for the first and last periods (i.e., i = 1 or M), the prior probability density function of the calendar ages of the early and late boundaries, say αi and βi, respectively, can be written as:
f(αi,βi|A,B)=1(A−B)2I(αi,βi|A,B), | (5) |
where I(αi,βi|A,B) is an indicator function describing the temporal order of the calendar ages of the early and late boundaries. It can be simply defined as:
I(αi,βi|A,B)={1,A>αi>βi>B0,otherwise. | (6) |
However, for the internal periods, the prior probability density function of the calendar ages of the early and late boundaries should be treated differently according to the relationship with the neighboring periods. Without loss of generality, for the ith internal period, where i=2,⋯,M−1, the prior probability density function of the early and late boundaries, say αi and βi, respectively, can be written as:
f(αi,βi|α1,βM)=1(α1−βM)2I(αi,βi|αi−1,βi−1,α1,βM), | (7) |
where I(αi,βi|αi−1,βi−1,α1,βM) is an indicator function describing the relationship of two neighboring periods. This mathematical framework is essentially the same as what Nicholls and Jones [14] proposed. Here three cases are considered: (1) "overlapping", which means the imbrication of two neighboring periods. An example is the initial coexistence of two archaeological cultures and subsequent replacement of one culture by another. Therefore, the indicator function can be expressed as:
I(αi,βi|αi−1,βi−1,α1,βM)={1,α1>αi−1>αi>βi−1>βi>βM0,otherwise, | (8) |
(2) "contiguous", which means a sudden transition from one period to another. An example is the rapid change of the depositional environment or the archaeological culture. As such, the indicator function can be defined as:
I(αi,βi|αi−1,βi−1,α1,βM)={1,α1>αi−1>βi−1=αi>βi>βM0,otherwise, | (9) |
and (3) "disjoint", which means that there is an interruption between two periods. An example is the existence of a depositional hiatus between two stratigraphical units or a gap of human occupation between two archaeological periods. As such, the indicator function can be written as:
I(αi,βi|αi−1,βi−1,α1,βM)={1,α1>αi−1>βi−1>αi>βi>βM0,otherwise, | (10) |
where " > " denotes "older than" and "=" denotes "same as".
The normal distribution is usually used for modeling 14C ages from the calibration curve [20,21]. However, the traditional normal distribution is less robust to outlying data. To better accommodate outliers in the 14C ages, we adopt the variance-multiplier normal distribution [22]. Let µ and ϵ be the mean and one standard deviation of a radiocarbon calibration curve, respectively. For the jth calendar age in the ith period, say θi,j, where i=1,2,⋯,M and j=1,2,⋯,Ni, the corresponding true 14C age is assumed to be normally distributed around the calibration curve mean, μμ(θi,j), confined by an additive variance of the laboratory 14C age and the calibration curve multiplied by a positive number. Hence, the true 14C age can be modeled as:
ˆri,j∼N(μ(θi,j),δω2i,j), | (11) |
where ω2i,j=σ2i,j+ϵϵ2(θi,j) is the sum of the laboratory 14C age variance (i.e., σ2i,j) and the calibration curve variance (i.e., ϵϵ2(θi,j)) and δ is a positive number. Let δ be a random variable with an improper prior probability density function defined as:
f(δ)=1δ, | (12) |
where δ∈(0,+∞). As such, for the ith period, where i=1,2,⋯,M, the likelihood function can be expressed as:
L(ri|θθi,σσi,δ)=∏Nij=11√2πδωi,jexp{−[ri,j−μμ(θi,j)]22δω2i,j}. | (13) |
An integrated likelihood can be obtained by integrating out hyperparameter δ (Supplementary), yielding:
L(ri|θθi,σσi)=Γ(Ni2)(∑Nij=1[ri,j−μμ(θi,j)]22ω2i,j)−Ni2, | (14) |
where Γ (·) denotes the gamma function.
Making use of Bayes' theorem and the integrated likelihood function (Eq. 14), for the ith period, where i=1,2,⋯,M, the joint posterior probability density function of parameters θi, αi, and βi, can be expressed as:
p(θθi,αi,βi|ri,σσi,A,B)∝L(ri|θθi,σσi)×f(θθi|αi,βi)×f(αi,βi|α1,βM)×f(α1,βM|A,B) |
=Γ(Ni2)(∑Nij=1[ri,j−μμ(θi,j)]22ω2i,j)−Ni2×1(αi−βi)NiI(θθi|αi,βi)×1(α1−βM)2I(αi,βi|αi−1,βi−1,α1,βM)×1(A−B)2I(α1,βM|A,B). | (15) |
Note that the prior indicator functions of the calendar ages and the boundaries have different expressions according to the ordering of the ages and the relationship of the boundaries as given above.
The structure of the above equation (Eq. 15) is so complex that we are unable to derive an analytical expression of the marginal posterior probability distribution function for each parameter. Therefore, this model framework can be implemented using the MCMC method. The Metropolis-approximated Gibbs sampler is used to generate a Markov Chain for the parameters one at a time, and a Metropolis-Hastings step is used to update each parameter [23]. Specifically, for the ith period of the sequence, where i=1,2,⋯,M, we choose to update the parameters in the order θi → αi → βi. Given the current state of the calendar ages of the 14C ages, say θi, and the early and late boundaries, say αi and βi, respectively, we first update θi using a uniform random walk process on the interval [αi, βi] such that:
θθ∗i|αi,βi=θθi+ρ×u, | (16) |
where u denotes a vector of Ni independent random numbers drawn from a uniform distribution on [–1, 1], and ρ is the step size of the random walk. Then, we update the calendar age of the early boundary of the ith period of the sequence, say αi, where i=1,2,⋯,M, using the triangular distribution. For the first period, the calendar age of the early boundary, say α1, can be updated according to:
α∗1|α1,A,θθ∗1∼T(max(θθ∗1),α1,A), | (17) |
where T (·) denotes the triangular distribution with three parameters describing the lower bound, the peak, and the upper bound, respectively. While for other periods (i.e., i > 1), special treatments are required, depending on the relationship of two neighboring periods. Three cases are considered accordingly. For the "overlapping" boundaries, the calendar age of the early boundary can be updated according to:
α∗i|αi,θθ∗i−1,θθ∗i∼T(max(θθ∗i),αi,min(θθ∗i−1(θθ∗i−1>αi))), | (18) |
while for the "contiguous" boundaries, the calendar age of the early boundary can be updated according to:
α∗i|αi,θθ∗i−1,θθ∗i∼T(max(θθ∗i),αi,min(θθ∗i−1)), | (19) |
and for the "disjoint" boundary, the calendar age of the early boundary can be updated according to:
α∗i|αi,βi−1,θθ∗i∼T(max(θθ∗i),αi,βi−1), | (20) |
Finally, we update the calendar age of the late boundary of the ith period, say βi, where i=1,2,⋯,M, using the triangular distribution too. Similarly, for the first M − 1 periods, three cases are considered. For the "overlapping" boundaries, the calendar age of the late boundary can be updated according to:
β∗i|βi,θθ∗i,θθ∗i+1∼T(max(θθ∗i+1(θθ∗i+1<βi)),βi,min(θθ∗i)); | (21) |
while for the "contiguous" boundaries, the calendar age of the late boundary, say βi, can be updated by simply setting:
β∗i=α∗i+1; | (22) |
and for the "disjoint" boundaries, the calendar age of the late boundary can be updated according to:
β∗i|βi,α∗i+1,θθ∗i∼T(α∗i+1,βi,min(θθ∗i)). | (23) |
For the last period (i.e., i = M), the calendar of the late boundary, say βM, can be updated according to:
β∗M|βM,B,θθ∗M∼T(B,βM,min(θθ∗M)). | (24) |
The acceptance or rejection of the proposed move of the chains is determined using the Metropolis-Hastings algorithm [24]. The state-of-the-art IntCal20 [25], SHCal20 [26], and Marine20 [27] calibration datasets are also provided for modeling of the 14C ages. For the sake of computing, all of the laboratory and the calibration curve radiocarbon ages are converted to the F14C space during the modeling process [20]. The convergence of the chains is monitored using the method proposed by Gelman [28]. Several chains are run parallelly with different starting points for each set of simulations, and the ratio of the between-sequence to with-sequence variance in terms of the potential scale reduction factor, ˆR is calculated. If ˆR is below the critical value of 1.1, the chains can be regarded as convergent. Once the chains for a model parameter converge, they are mixed to allow the estimate of the empirical probability density function that mimics the posterior probability distribution of this parameter. We implemented the above procedure in the Matlab® environment, yielding a software package named as MatCalib v2.0. The code is available on Mendeley Data [29] and also accessible in the Supplementary.
Running the software package is quite straightforward. A tab delimited spread sheet template is provided, which allow users to enter their data, including laboratory codes, depths, 14C ages and errors, reservoir ages and errors if any, specify the calibration curve (e.g., IntCal13 or IntCal20) for each age, and define the temporal ordering of the ages (e.g., "ordered", "unordered", or "coeval") in each period. The only file users need to modify is MatCalib_main.m, which provides a textual interface enabling users to read their data, define their project-specific problem (e.g., boundary relationships), specify parameters for MCMC simulations, and run the model to obtain, analyze, and visualize the results. In the Matlab environment, users need to open the main program MatCalib_main, give the nearest years that the modeled calendar ages are to be rounded to, choose the time scale in which the modeled calendar ages are to be reported, and specify the parameters for Monte Carlo simulations such as the number of chains to be run for convergence diagnoses as well as the length, burn-in period, and thinning interval of the chains. Then, MCMC simulations are conducted using the function age_modeling and the convergence of the chains is tested using the function convergence.
Upon completion of the MCMC simulations, the empirical probability density function, as well as the 68.2% and 95.4% HPD regions of the modeled calendar ages of the 14C ages and the early and late boundaries of each period, are estimated using the method of Lougheed and Obrochta [30], and the results are saved automatically to files for subsequent analyses. Users can plot the results either against sample numbers using the function plot_ages, or against depths using the function plot_age_depth if depth information is provided. Users also can compare any two ages by calculating their posterior difference and calculate the pooled mean of several ages using the functions age_difference and pooled_mean, respectively. The results are also given in terms of the empirical probability density function as well as the 68.2% and 95.4% HPD regions, and they can be plotted using the functions plot_difference and plot_pooled_mean, respectively.
It is noteworthy that the MCMC method is simulation-based, which requires a large number of iterations to produce reliable (i.e., converged and well mixed) results. Therefore, the model may yield slightly different results at each time of a run. Another caveat is that models may be so inconsistent with the data that no possible solutions will be yielded, particularly in the case of the presence of extreme outliers. Therefore, users should choose and remove the outliers before running the model to produce valid results. In other cases, the solution space may be very fragmentary, providing that the MCMC solutions will be different with each different starting position. Keeping these properties of the MCMC method in mind, users must check the reproducibility by making multiple runs of the model with a different number of iterations, burn-in period, thinning size, and initial value (automatically reset for each run using the function initialization).
Some of the functionalities of this software package are demonstrated using two datasets of 14C ages. One is obtained from the Neolithic archaeological sites on the Shandong Peninsula in North China, which is used to determine the timing of the Longshan and Yueshi cultures in this area. The other is retrieved from a sedimentary core on the southeastern Swedish Baltic coast, which is used to establish an age-depth relationship and to quantify the sedimentary hiatus induced by an erosional event during the transition from a lacustrine to a lagoonal environment. In both cases, a set of simulations with three chains are run for 20,000 iterations for each parameter. Every 20th iteration is kept after a burn-in period of 1000 iterations to reduce autocorrelation.
The dataset used in this study consists of nine 14C ages from sites belonging to the Longshan cultural period and 10 ages belonging to the Yueshi cultural period [31]. These ages are put together into a chronological sequence with two periods corresponding to the Longshan (Period 1) and Yueshi (Period 2), respectively. The latest Longshan site was dated to 3570 ± 80 14C yr BP and the earliest Yueshi site was dated to 3760 ± 145 14C yr BP, implying that these two cultures may have coexisted for a while [32]. Therefore, overlapping boundaries are used in the age modeling. As the 14C ages in each cultural period come from individual sites and their stratigraphical relationship is unknown, the temporal order of the ages is set to "unordered". The probability distributions of the modeled calendar ages are presented in Figure 3. The 95.4% HPD regions of calendar age of the boundaries are listed in Table 1 for the sake of comparison with those obtained from OxCal.
Software | Longshan Period (95.4% HPD region) | Yueshi Period (95.4% HPD region) | |||
Start | End | Start | End | ||
OxCal | 3040–2344 | 2332–1602 | 2552–1891 | 1693–1159 | |
MatCalib | 2885–2750 | 2005–1870 | 2490–2335 | 1480–1300 |
The Longshan culture, colloquially referred to as the black pottery culture, is a late-Neolithic culture in the middle and lower Yellow River areas of northern China from about 3000 to 1900 BCE. The culture flourished dramatically and expanded to the Shandong Peninsula at about 2885–2750 BCE as defined by the 95.4% HPD region (Table 1). It then decreased and vanished around 2005–1870 BC (95.4% HPD region). The Yueshi culture rose locally at about 2490–2335 BC (95.4% HPD region). It then replaced the Longshan culture and gradually evolved into the Bronze Age at about 1480–1300 BC (95.4% HPD region). These timings lie in the 95.4% age ranges obtained from OxCal (Table 1), suggesting that MatCalib can yield age estimates much firmer than OxCal. As defined by the age difference between the early boundary of Period 2 (i.e., Yueshi period) and the late boundary of Period 1 (i.e., Longshan period) (Figure 3), these two cultures may have coexisted for about 350–590 years (95.4% HPD region) on the Shandong Peninsula (Figure 4). Again, OxCal yields a relatively large overlap (280–740 years, 95% HPD region) of these two cultures.
Previous work about the 14C dating of a sediment core from the Hunnemara ancient lake (14°53' E, 56°10' N) on the southeastern Swedish Baltic coast revealed a depositional hiatus across the transition from the lacustrine to the lagoonal facies [33], which may have resulted from erosion during the rapid sea-level rise about 7600 years ago [34]. The lake is developed within a closed basin with an outflow threshold around 3.0 m as1. It has been drained and converted to agricultural land and a meadow by local farmers in the 19th century and is now partly used for a garbage dump. An 8.5-m-long core was taken near the center of the lake basin. A total of 10 14C ages on bulk sediments were obtained, which are incorporated into a chronological sequence consisting of two periods representing the lacustrine (Period 1) and lagoonal (Period 2) environments, respectively. Period 1 comprises three 14C ages, whereas Period 2 contains seven 14C ages subject to a local reservoir age offset of −108 ± 24 years. Therefore, the Marine20 calibration curve is applied on these ages [27]. As the ages are obtained from different depths along the core, the law of superposition was strictly enforced by using a temporal order constraint (e.g., "ordered") on the ages in both periods. The hiatus lies at 575 cm, which is marked by a non-depositional surface separating the two periods. Therefore, a "disjoint" relationship of the internal boundaries is applied.
The modeled calendar ages of the radiocarbon ages are plotted against their depths (Figure 5). The age-depth relationship reveals variable sedimentation rates within the lake basin. The results suggest that the lacustrine environment (Period 1) prevailed from 11,745–11,545 BP (95.4% HPD region) to 8530–8145 BP (95.4% HPD region), while the lagoonal environment (Period 2) occurred in the basin at about 7410–7040 BP (95.4% HPD region) as a result of rapid sea-level rise [34] and terminated at about 3090–2925 BP (95.4% HPD region). As defined by the posterior age difference of the late boundary of Period 1 and the early boundary of Period 2, the sedimentary hiatus is estimated to be 870–1355 years (95.4% HPD region) (Figure 6).
Bayesian methods are widely used for modeling 14C ages. It represents a flexible yet robust approach to integrating prior knowledge about the ages and their temporal or stratigraphical context, thereby yielding a more precise timing of the geological or archaeological events than the traditional methods [8]. Although some key components of the Bayesian methods for 14C age modeling have already been well explained [8,14,21], a formal treatment of several other outstanding issues (e.g., the relationship of neighboring boundaries) is lacking. Here a complete and detailed description of the analytical framework is presented along with a Matlab package for software implementation. Therefore, it can be used not only as an introductory guild to, but also a tutorial on Bayesian 14C age modeling.
There exist several software packages designed to perform Bayesian 14C age modeling [11,17,18]. Most of which shares some common features. For example, the law of superposition is the backbone of all Bayesian methods, which has been strictly enforced to observe the constraint of temporal order. Yet the model presented here differs from these tools in several aspects. First is the treatment of the relationship between two neighboring periods. As illustrated by the first example, coexistence of two cultures appears to be a common phenomenon in archaeology [32]. Similar to OxCal, an overlapping relationship is considered in this model, while this functionality was not implemented in other tools. Second is the application of different calibration curves on 14C ages of different types. This is particularly necessary if ages were obtained from different depositional environments (e.g., marine versus terrestrial). As illustrated by the second example, ages in the first period are obtained from a terrestrial setting, whereas those in the second period are retrieved from a marine environment subject to the reservoir effect [35]. As such, a local reservoir age and the marine calibration curve should be applied on these ages. Other minor differences of this model from the existing tools include the use of the non-Gaussian likelihood function to accommodate outlying data, the implementation of convergence diagnosis of the Markov Chains, and the employment of a triangular distribution to update the boundaries during the simulations.
The accompanying Matlab package (i.e., MatCalib 2.0) contains a wide range of functionality for Bayesian statistical analyses of 14C ages, including the inference of the calendar age of the individual 14C ages and the boundary of the periods. Random samples of these parameters generated from the MCMC simulations are saved automatically to a data file for future use without running the model again. Based on these data, the empirical probability density functions of these parameters are built up first and then the HPD regions and the median-probability age are calculated. Further analyses such as pooled mean of ages, timing, and duration of an event can be conducted and visualized using the high-quality graphical output (e.g., Figures 3–6). In future, the software will be improved to enable the modeling of multi-source ages, for example, calendar ages with errors and/or tie points for events without errors. For these age types, a calibration curve is no longer needed. If depth information is available, an age-depth relationship may be established through polynomial regression (e.g., CLAM [36]), Brownian bridge (e.g., OxCal [37]), or interpolation (e.g., Undatable [38]) on the stored random samples at known depths, instead of employing a stochastic depositional model such as Bchron [39] and Bacon [40].
The traditional method of 14C age calibration maps individual 14C ages to the calibration curve to obtain the calendar age. However, owing to the large fluctuations of the atmospheric 14C level, the resulting calibrated age is not a point estimate; rather, it is usually treated as a random variable expressed in terms of a probability distribution with one or several HPD regions. Most importantly, the probability distribution of the calibrated age could be extremely broad if the radiocarbon age hits a radiocarbon plateau in the calibration curve, thereby yielding several large HPD regions. The large uncertainty in the calibrated age hinders the precise timing of geological and/or archaeological events. Conceptually different from the traditional method, Bayesian model provides such an analytical framework that integrate prior knowledge about the ages and their temporal order as a constraint, which in turn can significantly narrow down the uncertainty of the calibrated ages.
The hierarchical Bayesian model presented in this paper provides a flexible approach to the precise timing of archaeological and/or geological events by making use of a suite of 14C ages and their stratigraphical context (temporal order). Implemented in the Matlab environment, the accompanying software package is not only an extension to the existing MatCal package that is only able to calibrate individual 14C ages [30], but also an alternative or complement to the online calibration tools such as BCal [18] and OxCal [11]. Also, this model can be used as a building block for age-depth modeling. If the depth information of the 14C ages is available, a probabilistic expression of the age-depth relationship can be readily obtained through either bootstrap regression or Brownian bridge transformation.
This work was supported by the National Science Foundation of China (grant number 41971102). My gratitude is due to the three anonymous reviewers for their insightful comments and constructive suggestions that greatly improved the manuscript.
The author declares no conflict of interest.
1. Derivation of the integrated likelihood function.
2. Matlab code for implementing Bayesian radiocarbon age modeling.
[1] |
Bartel DP (2004) MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116: 281-297. doi: 10.1016/S0092-8674(04)00045-5
![]() |
[2] |
Kim YJ, Zheng B, Yu Y, et al. (2011) The role of mediator in small and long noncoding RNA production in Arabidopsis thaliana. EMBO J 30: 814-822. doi: 10.1038/emboj.2011.3
![]() |
[3] |
Zhang BH, Pan XP, Wang QL, et al. (2005) Identification and characterization of new plant microRNAs using EST analysis. Cell Res 15: 336-360. doi: 10.1038/sj.cr.7290302
![]() |
[4] |
Hammond SC, Bernstein E, Beach D, et al. (2000) An RNA-directed nuclease mediates posttranscriptional gene silencing in Drosophila cells. Nature 404: 293-296. doi: 10.1038/35005107
![]() |
[5] |
Kidner CA, Martienssen RA (2005) The developmental role of microRNA in plants. Curr Opin Plant Biol 8: 38-44. doi: 10.1016/j.pbi.2004.11.008
![]() |
[6] |
Baloch IA, Barozai MYK, Din M (2013) MicroRNAs: the mega regulators in eukaryotic genomes. Pure Appl Biol 2: 83-88. doi: 10.19045/bspab.2013.23002
![]() |
[7] |
Bai M, Yang GS, Chen WT, et al. (2012) Genome-wide identification of Dicer-like, Argonaute and RNA dependent RNA polymerase gene families and their expression analyses in response to viral infection and abiotic stresses in Solanum lycopersicum. Gene 501: 52-62. doi: 10.1016/j.gene.2012.02.009
![]() |
[8] | Barozai MYK (2012) Insilico identification of microRNAs and their targets in fiber and oil producing plant Flax (Linum usitatissimum L.). Pak J Bot 44: 1357-1362. |
[9] |
Gao P, Bai X, Yang L, et al. (2011) Osa-MIR393: a salinity- and alkaline stress-related microRNA gene. Mol Biol Rep 38: 237-242. doi: 10.1007/s11033-010-0100-8
![]() |
[10] |
Shui XR, Chen ZW, Li JX (2013) MicroRNA prediction and its function in regulating drought-related genes in cowpea. Plant Sci 210: 25-35. doi: 10.1016/j.plantsci.2013.05.002
![]() |
[11] |
Xie FL, Huang SQ, Guo K, et al. (2007) Computational identification of novel microRNAs and targets in Brassica napus. FEBS Lett 581: 1464-1474. doi: 10.1016/j.febslet.2007.02.074
![]() |
[12] |
Zhang BH, Pan XP, Stellwag EJ (2008) Identification of soybean microRNAs and their targets. Planta 229: 161-182. doi: 10.1007/s00425-008-0818-x
![]() |
[13] |
Barozai MYK, Irfan M, Yousaf R, et al. (2008) Identification of micro-RNAs in cotton. Plant Physiol Biochem 46: 739-751. doi: 10.1016/j.plaphy.2008.05.009
![]() |
[14] |
Zhang BH, Wang QL, Wang KB, et al. (2007) Identification of cotton microRNAs and their targets. Gene 397: 26-37. doi: 10.1016/j.gene.2007.03.020
![]() |
[15] |
Zhang B, Pan X, Cannon CH, et al. (2006) Conservation and divergence of plant microRNA genes. Plant J 46: 243-259. doi: 10.1111/j.1365-313X.2006.02697.x
![]() |
[16] |
Frazier TP, Xie F, Freistaedter A, et al. (2010) Identification and characterization of microRNAs and their target genes in tobacco (Nicotiana tabacum). Planta 232: 1289-1308. doi: 10.1007/s00425-010-1255-1
![]() |
[17] |
Xie F, Frazier T, Zhang B (2010) Identification and characterization of microRNAs and their targets in the bioenergy plant switchgrass (Panicum virgatum). Planta 232: 417-434. doi: 10.1007/s00425-010-1182-1
![]() |
[18] |
Barozai, MYK, Din M, Baloch IA (2013) Structural and functional based identification of the bean (Phaseolus) microRNAs and their targets from Expressed Sequence Tags. J Struct Funct Genomics 14: 11-18. doi: 10.1007/s10969-013-9152-z
![]() |
[19] |
Din M, Barozai MYK (2014) Profiling microRNAs and their targets in an important fleshy fruit: Tomato (Solanum lycopersicum). Gene 535: 198-203. doi: 10.1016/j.gene.2013.11.034
![]() |
[20] | Din M, Barozai MYK (2014) Profiling and characterization of eggplant (Solanum melongena L.) microRNAs and their targets. Mol Biol Rep 41: 889-894. |
[21] | Din M, Barozai MYK, Baloch IA (2016) Profiling and annotation of microRNAs and their putative target genes in chilli (Capsicum annuum L.) using ESTs. Gene Rep 5: 62-69. |
[22] |
Muchero W, Diop NN, Bhatetal PR (2009) A consensus genetic map of cowpea (Vigna unguiculata (L) Walp) and synteny based on EST-derived SNPs. Proc Natl Acad Sci U.S.A. 106:18159-18164. doi: 10.1073/pnas.0905886106
![]() |
[23] |
Pule-Meulenberg F, Belane AK, Krasova-Wade T, et al. (2010) Symbiotic functioning and bradyrhizobial biodiversity of cowpea (Vigna unguiculata L. Walp) in Africa. BMC Microbiol 10: 89. doi: 10.1186/1471-2180-10-89
![]() |
[24] | Griffiths-Jones S (2004) The microRNA registry. Nucleic Acids Res 32D: 109-111. |
[25] |
Altschul SF, Gish W, Miller W, et al. (1990) Basic local alignment search tool. J Mol Biol 215: 403-410. doi: 10.1016/S0022-2836(05)80360-2
![]() |
[26] |
Altschul SF, Madden TL, Schäffer AA, et al. (1997) Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25: 3389-3402. doi: 10.1093/nar/25.17.3389
![]() |
[27] |
Zuker M (2003) Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res 31: 3406-3415. doi: 10.1093/nar/gkg595
![]() |
[28] |
Barozai MYK (2012) Identification and characterization of the microRNAs and their targets in Salmo salar. Gene 499: 163-168. doi: 10.1016/j.gene.2012.03.006
![]() |
[29] | Ambros V, Bartel B, Bartel DP, et al. (2003) A uniform system for microRNA annotation. RNA 9: 277-279. |
[30] |
Dai X, Zhao PX (2011) psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res 39: 155-159. doi: 10.1093/nar/gkq766
![]() |
[31] |
Barozai MYK (2012) The microRNAs and their targets in the channel catfish (Ictalurus punctatus). Mol Biol Rep 39: 8867-8872. doi: 10.1007/s11033-012-1753-2
![]() |
[32] |
Kruger J, Rehmsmeier M (2006) RNAhybrid: microRNA target prediction easy, fast and flexible. Nucl Acids Res 34: 451-454. doi: 10.1093/nar/gkj455
![]() |
[33] | Barozai MYK, Husnain T (2011) Identification of biotic and abiotic stress up-regulated ESTs in Gossypium arboreum. Mol Biol Rep 39: 1011-1018. |
[34] | Barozai MYK, Wahid AH (2012) In silico identification and characterization of cumulative abiotic stress responding genes in Potato (Solanum tuberosum L.). Pak J Bot 44: 57-69. |
[35] |
Barozai MYK, Kakar AG, Din M (2012) The relationship between codon usage bias and salt resistant genes in Arabidopsis thaliana and Oryza sativa. Pure Appl Biol 1: 48-51. doi: 10.19045/bspab.2012.12005
![]() |
[36] | Barozai MYK, Kakar S, Sarangzai AM (2013) Profiling the carrot (Daucus carota L.) microRNAs and their targets. Pak J Bot 45: 353-358. |
[37] |
Wang J, Yang X, Xu H, et al. (2012) Identification and characterization of microRNAs and their target genes in Brassica oleracea. Gene 505: 300-308. doi: 10.1016/j.gene.2012.06.002
![]() |
[38] | Barozai MYK (2013) Identification of microRNAs and their targets in Artemisia annua L. Pak J Bot 45: 461-465. |
[39] |
Ghani A, Din M, Baloch IA, et al. (2013) Identification of MicroRNA in 12 plant species of fabaceae. Pure Appl Bio 2: 104-115. doi: 10.19045/bspab.2013.23005
![]() |
[40] | Orlov YL, Dobrovolskaya O, Yuan CH, et al. (2012). Integrative computer analysis of antisense transcripts and miRNA targets in plant genomes. J Stress Physiol Biochem 8: S7. |
[41] |
Barozai MYK (2012) The novel 172 sheep (Ovis aries) microRNAs and their targets. Mol Biol Rep 39: 6259-6266. doi: 10.1007/s11033-012-1446-x
![]() |
[42] |
Chen L, Ren YY, Zhang YY, et al. (2012) Genome-wide identification and expression analysis of heat-responsive and novel microRNAs in Populus tomentosa. Gene 504: 160-165. doi: 10.1016/j.gene.2012.05.034
![]() |
[43] |
Ji Z, Wang G, Xie Z, et al. (2012) Identification and characterization of microRNA in the dairy goat (Capra hircus) mammary gland by Solexa deep sequencing technology. Mol Biol Rep 39: 9361-9371. doi: 10.1007/s11033-012-1779-5
![]() |
[44] |
Barozai MYK (2012) The microRNAs and their targets in the channel catfish (Ictalurus punctatus). Mol Biol Rep 39: 8867-8872. doi: 10.1007/s11033-012-1753-2
![]() |
[45] |
Yu J, Wang F, Yang GH, et al. (2006). Human microRNA clusters: genomic organization and expression profile in leukemia cell lines. Biochem Biophys Res Commun 349: 59-68. doi: 10.1016/j.bbrc.2006.07.207
![]() |
[46] |
Jones-Rhoades MW, Bartel DP (2004) Computational identification of plant microRNAs and their targets, including a stress induced miRNA. Mol Cell 14: 787-799. doi: 10.1016/j.molcel.2004.05.027
![]() |
[47] |
Crooks GE, Hon G, Chandonia JM, et al. (2004) Web-Logo: a sequence logo generator. Genome Res 14:1188-1190. doi: 10.1101/gr.849004
![]() |
[48] |
Larkin MA, Blackshields G, Brown NP, et al. (2007) ClustalW and ClustalX version 2. Bioinform 23: 2947-2948. doi: 10.1093/bioinformatics/btm404
![]() |
[49] | Zeng CY, Wang WQ, Zheng Y, et al.(2009) Conservation and divergence of microRNAs and their functions in Euphorbiaceous plants. Nucleic Acids Res 38: 981-995. |
[50] |
Bartel DP (2009) MicroRNAs: target recognition and regulatory functions. Cell 136: 215-233. doi: 10.1016/j.cell.2009.01.002
![]() |
[51] | Kohli P, Kalia M, Gupta R (2015) Pectin Methylesterases: A Review. J Bioprocess Biotech 5: 228. |
[52] |
Whitney SM, Andrews TJ (2001) The gene for the ribulose-1, 5-bisphosphate carboxylase/oxygenase (Rubisco) small subunit relocated to the plastid genome of tobacco directs the synthesis of small subunits that assemble into Rubisco. Plant Cell 13: 193-205. doi: 10.1105/tpc.13.1.193
![]() |
[53] |
Ballester AR, Molthoff J, de Vos R, et al. (2010) Biochemical and molecular analysis of pink tomatoes: deregulated expression of the gene encoding transcription factor SlMYB12 leads to pink tomato fruit color. Plant Physiol 152: 71-84. doi: 10.1104/pp.109.147322
![]() |
[54] |
Kodaira KS, Qin F, Tran LS, et al. (2011) Arabidopsis Cys2/His2 zinc-finger proteins AZF1 and AZF2 negatively regulate abscisic acid-repressive and auxin-inducible genes under abiotic stress conditions. Plant Physiol 157: 742-756. doi: 10.1104/pp.111.182683
![]() |
[55] |
Soria-Guerra RE, Rosales-Mendoza S, Gasic K, et al. (2011) Gene expression is highly regulated in early developing fruit of apple. Plant Mol Biol Rep 29: 885-897. doi: 10.1007/s11105-011-0300-y
![]() |
[56] |
Yadav SK (2010) Cold stress tolerance mechanisms in plants. A review. Agron Sustain Dev 30: 515-527. doi: 10.1051/agro/2009050
![]() |
[57] | Prasad PVV, Staggenborg SA, (2008) Impacts of drought and/or heat stress on physiological, developmental, growth, and yield processes of crop plants, In Ristic, Z. Author, Response of Crops to Limited Water. Madison, WI, USA, 301-355. |
[58] | Qados AMSA (2011) Effect of salt stress on plant growth and metabolism of bean plant Viciafaba (L.). J Saudi Soc Agric Sci 10: 7-15. |
[59] |
Rejeb IB, Pastor V, Mauch-Mani B (2014) Plant Responses to Simultaneous Biotic and Abiotic Stress: Molecular Mechanisms. Plants 3:458-475. doi: 10.3390/plants3040458
![]() |
[60] | Sheshadri SA, Nishanth MJ, Simon B (2016) Stress-mediated cis-element transcription factor interactions interconnecting primary and specialized metabolism in planta. Front Plant Sci 7: 1725. |
[61] | Fluhr R (2001) Sentinels of disease. Plant resistance genes. Plant Physiol 127: 1367-1374. |
[62] |
Umezawa T, Yoshida R, Maruyama K, et al. (2004) SRK2C, a SNF1-related protein kinase 2, improves drought tolerance by controlling stress-responsive gene expression in Arabidopsis thaliana. Proc Natl Acad Sci U.S.A. 101: 17306-17311. doi: 10.1073/pnas.0407758101
![]() |
[63] |
Narayan A, Sachdeva P, Sharma K, et al. (2007) Serine threonine protein kinases of mycobacterial genus: phylogeny to function. Physiol genomics 29: 66-75. doi: 10.1152/physiolgenomics.00221.2006
![]() |
Software | Longshan Period (95.4% HPD region) | Yueshi Period (95.4% HPD region) | |||
Start | End | Start | End | ||
OxCal | 3040–2344 | 2332–1602 | 2552–1891 | 1693–1159 | |
MatCalib | 2885–2750 | 2005–1870 | 2490–2335 | 1480–1300 |