A note on the population genetic consequences of delayed larval development in insects

EVOLUTIONARY GENETICS

RESEARCH ARTICLE

A note on the population genetic consequences of delayed larval development in insects

Marcos Mattoso de Salles; Paulo A. Otto

Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, SP, Brazil

^{Send correspondence to}

ABSTRACT

Observations by Dobzhansky's group in the 1940s suggesting that the presence of recessive genotypes could account for lower larval developmental rates in Drosophila melanogaster were not confirmed at the time and all subsequent investigations on this subject focused on the analysis of ecological models based on competition among pre-adult individuals. However, a paper published in this journal in 1991 eventually confirmed the finding made by Dobzhansky and his co-workers. In this report, we provide a theoretical analysis of the population genetic effects of a delay in the rate of larval development produced by such a genetic mechanism.

Keywords: difference and differential equations, Drosophila, dynamic systems, larval development, population genetics.

Introduction

Based on an analysis of samples from natural populations of Drosophila pseudoobscura collected from three localities on Mount San Jacinto (California), Dobzhansky et al. (1942) showed that homozygotes for genes located on the second and third chromosomes in this species had a lower developmental rate than expected. These authors suggested that, under controlled culture conditions, the viability and developmental rates were dependent on the demographic density of the larvae. In contrast to the genetic finding, the latter observation was subsequently confirmed by several authors (Bakker, 1959, 1969; De Witt, 1960; Robertson, 1963, 1964; Barker and Podger, 1970; Huang et al., 1971; Tosic and Ayala, 1981; Mather and Caligari, 1981; Ménsua and Moya, 1983; Adell et al., 1988), not only for drosophilid populations but for other insects as well. As a result, all subsequent research on this subject has centered on the analysis of ecological models based on competition among pre-adult individuals. Oliveira and Cordeiro (1981) detected a major gene effect on delayed development in Drosophila melanogaster that was later associated with a recessive gene located on the second chromosome of this fly (Oliveira et al., 1991).

Since the development of molecular genetics in the early 1990s many authors have examined the role of genes involved in the regulation of developmental timing and body size determination in Drosophila. Most of these studies focused on genes related to apoptosis and to the production or release of ecdysone and other ecdysteroids (molting hormones). The results of some relevant articles (all dealing with this hormone in D. melanogaster) are summarized below. Other authors described the importance of additional gene-controlled systems, such as those involved in the regulation of adenosine deaminase-related growth factor ADGF-A, in developmental timing (Dolezal et al., 2005) and in the regulation of the activities of oxidative phosphorylation complexes through nuclear-encoded mitochondrial tyrosyl-tRNA synthetases (Meiklejohn et al., 2013). Oldham et al. (2000) and Quinn et al. (2012) provide appropriate reviews of this subject.

Sliter and Gilbert (1992) showed that loss of function mutations at the dre4 gene caused ecdysteroid deficiency and developmental delay or arrest. Later, Bialecki et al. (2002) showed that isoform-specific null mutants for the orphan member E75A gene had a reduced ecdysteroid titer that resulted in developmental delay. McBrayer et al. (2007) showed that prothoracicotropic (PTTH) hormone production was not essential for molting or metamorphosis but that its lack resulted in significantly delayed larval development. Ghosh et al. (2010) showed that the gap gene giant (gt) regulated ecdysone production through specification of PTTH and that mutants for this gene exhibited delayed larval development (in addition to other effects).

The purpose of the present theoretical investigation was to analyze the population genetic effects of a delay in the larval developmental rate produced by an autosomal recessive mode of inheritance.

Using elements from the theory of finite differences, difference and differential equations, and dynamic systems, we examined the problem based on a model with discrete generations and with the total number of individuals (adults plus larvae) kept constant. We also summarize the results obtained with alternative models involving an increase or decrease in the total population number. The case of continuous generations (to be presented in detail elsewhere) is also briefly discussed.

Description of the main model in which the total number of individuals is kept constant over discrete generations

The central principle of this model is that the total number of individuals (adults + larvae) is kept constant in a population with discrete generations. For this, let A and a be a pair of segregating alleles at an autosomal locus. In the homozygous state the recessive gene a induces a delay in larval development (increase in the diapause phase), measured in an integer number of generations: crosses therefore take place randomly (panmixia) only among individuals belonging to the same generation. The population size, kept rigorously constant throughout generations, is assumed to be so large that the effects of random sampling genetic drift are negligible. Consequently, the mathematical treatment is deterministic and gene frequencies are interpreted as probabilities. The effects of differential migrations, selection pressure and mutation are also considered negligible.

In the equations describing the model we will use the following parameters: N₀= number of adults in the initial population at generation 0; this number equals the total number of individuals (adults plus larvae) at any time or generation K > 0; q₀= initial frequency of the recessive allele a [q₀= f₀(a)]; p₀= initial frequency of the normal allele A[p₀ = f₀(A) = 1-q₀]; T = delay in larval development determined by the genotype aa, measured in an integer number of generations; K = index corresponding to an integer number of generations; N_K= number of adult individuals in the population at the K-th generation; N'_K= number of larvae in the population at the K-th generation; q_K= frequency of the recessive gene a among adults at the K-th generation; p_K= frequency of the normal allele, also among adults belonging to the K-th generation; Q_K= frequency of the recessive gene in the total population (adults plus larvae) at the K-th generation; P_K= frequency of the normal allele, also in the total population, at the K-th generation (P_K = 1Q_K).

Derivation of pertinent equations

Based on the foregoing definitions:

(a) the total number of adults at generation K is given by

N_K =N_K-1(1-q_K-1²) if K<T or

N_K = N_K-1(1-q_K-1²)+N_K-1-Tq_K-1-T² if K > T;

(b) the genotypic proportions among AA, Aa and aa adult individuals at generation K are, respectively:

(N_K-1/N_K)p_K-1²,(N_K-1/N_K)2p_K-1q_K-1, and 0 if K<T or

(N_K-1/N_K)p_K-1²,(N_K-1/N_K)2p_K-1q_K-1, and (N_K-1-T/N_K)q_K-1-T² if K > T;

(c) the number of larvae aa present at generation K after the maturation of normal individuals is given by

N'_K = N₀q₀² +N₁q₁² +...+N_K-1q_K-1² if K < T or

N'_K = N_K-Tq_K-T² + N_K-T+1q_K-T+1² +...+ N_K-1q_K-1² if K > T.

Using the foregoing equations, it is not difficult to show that N_K+1 + N'_K+1 = N_K +N'_K = .. .= N₀ + 0 = N₀ and that, if P_K = f(A) in the total population, then P_K = (N_K/N₀)p_K = (N_K/N₀)(N₀p₀/N_K) = p₀, i.e., allele frequencies and the number of individuals remain constant when the whole population (adults plus larvae) is considered, as assumed by the model.

For K < T, we have, among adults of the population, a mechanism similar to the classic model of total selection against homozygous recessive individuals aa, despite the absence of selection against any individual of the population (when adults and larvae are considered together), as shown by the equations derived in the previous paragraph. This results in q_K+1 = q_K/(1+q_K), a recurrence equation that has the general solution in simple analytical form q_K = q₀/(1+Kq₀).

For K > T, we have

N_K+1 = N_K(1-q_K²)+N_K-Tq_K-T² =

N_K(1-q_K)(1+q_K)+N_K-Tq_K-T² =

N_Kp_K(2-p_K)+N_K-T(1-p_K-T)².

Since p_K+1 = (N_Kp_K²+N_Kp_Kq_K)/N_K+1 = N_Kp_K/N_K+1and, hence, N_K = N₀p₀/p_K, it follows that N₀p₀/p_K+1= N₀p₀(2-p_K)+N₀p₀(1-p_K-T)²/p_K-Tand

1/p_K+1 = 1/p_K-T + p_K-T -p_K,

p_K+1/p_K+1 = p_K-T + 1/p_K-T,

or, since p_K =N₀p₀/N_K, we also obtain

NK₊₁ = N_K-T + (N₀p₀)²(1/N_K-T - 1/N_K).

Dynamic behavior of allele frequencies

Below, we analyze the behavior of allele frequencies p_Kand q_K.

For K < T, p_K =(p₀ +Kq₀)/(1 + Kq₀), i.e., p_Kis a monotonically increasing function of K, whereas for K > T, 1/p_K+1 = p_K-T + 1/p_K-T - p_K.

Analysis of the function f(x)=x+1/x in the interval 0 <x<1 shows that f'(x) =1-1/x² < 0; f(x) is therefore a monotonically decreasing function of x in this interval.

If p_K > p_K-T, then we have p_K + 1/p_K < p_K-T + 1/p_K-T = p_K + 1/p_K+1and, therefore, p_K+1 < p_K; similarly, p_K <p_K-Timplies that p_K+1 > p_K, and p_K = p_K-Timplies that p_K+1 = p_K.

In the interval 0 < K < T, pK increases such that when K surpasses T, we have p_K> p_K-T; p_Kthen decreases to a new point where p_K < p_K-Tis reached and we will then have another increase in p_K, and so on. Thus, there will be an oscillation in gene frequencies among adults of the population. With the passing of generations, i.e., when K increases, the variations become smaller. In the limit case, when K tends to infinity, the gene frequencies reach an equilibrium (p_inf =p_e,q_inf =q_e).

Figures 1-3 show numerical examples of oscillatory convergence to equilibrium points (p_e) over 10 generations of random crosses, for populations starting with initial gene frequencies p₀= 0.1, ... , 0.9 and T values of T=1, 2, and 3. These graphs and others shown in this paper were prepared using routines and packages from Mathematica^® v.8.0 (Wolfram Research).

The preliminary global analysis presented above already shows a point of real theoretical interest, namely, that despite intense selective pressure among adult individuals (but not among the total population of adults plus larvae), the recessive allele is never eliminated from the population, but rather is kept in equilibrium at a point q_e> 0. This is true except when the number of generations involved in the developmental delay is very large or tends to infinity such that allele a is practically eliminated from the adult population and all larvae fail to develop.

Equilibrium gene frequencies

In a generic instant we have the following population composition (Figure 4).

If we consider an equilibrium situation where

N_K+1 = N_K = N_K-1 = N_K-2 = ... = N_K-T = N_e

q_K+1 = q_K = q_K-1 = q_K-2 = ... = q_K-T = q_e = 1-p_e,

from N_K + N'_K = N₀it follows that N_e +TN_eq_e² = N₀.

Since N_e = N₀p₀/p_e, it follows that N₀p₀(1+Tq_e²)/(1-q_e)=N₀, and, hence, p₀ +Tp₀q_e² = 1-q_eand, finally,

q_e = [(1+4Tp₀q₀)^1/2-1]/(2Tp₀) =

2q₀/[(1+4Tp₀q₀)^1/2+1]

.

As expected, q_e = q₀if T=0 (simple case of panmixia without selection) and q_e = 0 if T = ∞ (special case of complete selection against recessive individuals aa).

Considering that:

p_K =N_K-1p_K-1/N_K =N₀p₀/N_Kfor any K,

N_K = N_K-1(1-q_K-1²) for any K<T,

and

N_K = N_K-1(1-q_K-1²) + N_K-T-1q_K-T-1² for any K > T,

it is not difficult to show that the recurrence equation 1/p_K+1 = p_K-T + 1/p_K-T - p_Kcan be rewritten in the more convenient form

1/p_K = 1/p₀ -Σ_[j=1,T]{(1-p_K-j)²/p_K-j}.

This recurrence equation allows straightforward derivation of the equilibrium point p_e. Indeed, when K tends to infinity, the limit for the right-hand side of the above equation is 1/p₀ -T(1-p_e)²/p_eand it follows that p_e =p₀(1+Tq_e²) and p_e =1-q_e = 1 -[(1+4Tp₀q₀)^1/2-1]/(2Tp₀).

Table 1 shows the sets of equilibrium points q_efor several values of T, as a function of the initial values q₀= 0.0, 0.1, ..., 0.9, 1.0.

Figure 5 shows the sets of equilibrium values q_eas a function of the initial frequency values q₀and T varying from 0 to 15.

For the case T = 1, the recurrence equation above can be rewritten as

1/p_K+1 = 1/p₀ -(1-p_K)²/p_K or p_K+1 =

p₀p_K/[p_K-p₀(1-p_K)²]

.

The eigenvalue r responsible for the convergence of the series {p₀, ..., p_e} is obtained by directly differentiating p_K+1 = f(p_K) at equilibrium point p_e:

r = dp_K+1/dp_K|p_K = p_e =

-p₀² (q_e² -2q_e)/(p_e-p₀q_e²)² = p_e² - 1.

Since |r| = |p_e²-1| < 1 for any value of p_e in the interval 0 < p_e < 1, the equilibrium point p_eis asymptotically stable and convergence to this point occurs at a geometric rate. Since r is always smaller than zero in this interval, convergence to equilibrium is oscillatory.

Figure 6 shows the values of r and q_eas functions of the initial gene frequency q₀for the case T=1.

For a generic value of T the method that follows is used for determining convergence rates as well as equilibrium properties.

Equilibrium dynamics

As shown above, for any K > T,

1/p_K = 1/p₀ - Σ_{[j = 1,T]}{(1-p_K-j)²/p_K-j} or

p_K = 1/[1/p0 - Σ_{[j = 1,T]}{(1-p_K-j)²/p_K-j}].

We then define a set of T variables as follows:

x₁(K) p_K

x₂(K) p_K-1

. . X_K=

( . )

= ( . )

.

.

.

.

x_T(K)

p_K-T+1

Since

p_K+1 = [1/p₀ - Σ _{[j = 1,T]}{(1-p_K-j+1)²/p_K-j+1}] - 1

we have:

x₁(K+1) =

f[x₁(K), x₂(K), x₃(K+1), x₄(K), ... , x_T(K)]

x₂(K+1) =

x₁(K)

x₃(K+1) =

x₂(K)

x₄(K+1) =

x₃(K)

...........................................................................................................

x_T(K+1) =

x_T-1(K).

Considering the equilibrium point, at which

x₁(K)

p_e

x₂(K) p_e

.

. X_e=

( . )

=

( . )

.

.

. .

x_T(K)

p_e

we have, in an infinitesimal neighborhood of this point, DX(K+1) = A.DX(K), or, in expanded form,

Dx₁

a a a ... a a Dx₁

Dx₂

1 0 0 ... 0 0 Dx₂

Dx₃

0 1 0 ... 0 0 Dx₃

(.)_K+1=

(............................................................................) (.)_K

.

............................................................................ .

.

............................................................................ . Dx_T 0 0 0 ... 1 0 Dx_T

a = df/dx₁(K) |x=x_e = df/dx₂(K) |x=x_e = ... =

= df/dx_j(K) |x=x_e =

=p_e² -1.

The eigenvalues r₁, r₂, ... , r_Tare the solutions of the characteristic equation

det(A-rI) = (-1)^T(1+r+r²+ ... +r^T-1 - r^T/a) = 0

or

r^T /a=1 + r + r² + r³ + ... + r^T-2 + r^T-1,

which can also be rewritten as r^T/a = (1-r^T)/(1-r), r ≠ 1.

Based on the function f(r)=r^T+1 - (a+1)r^T + a, r ≠ 1, -1 < a < 0, it can be shown that f(0) =a<0, f(1) = 0 and df(r)/dr = f'(r) = (T+1)r^T-T(a+1)r^T-1. Also, as r tends to +∞, f(r) tends to +inf and as r tends to -∞, f(r) tends to +∞ if T is odd, and to -∞ if T is even. When T is even, f'(r) is an increasing function of r in the open interval (-∞, 0), decreasing in the open interval (0, (a+1)T/(T+1)) and again increasing for r > 0. Therefore, r = 0 is a relative maximum and r = (a+1)T/(T+1) a relative minimum. When T is odd, f'(r) is a decreasing function of r in the open interval (-inf, (a+1)T/(T+1)) and an increasing function of r for r> (a+1)T/(T+1). Therefore, r=0 in the case when T is odd is an inflection point, whereas r = (a+1)T/(T+1) is a relative minimum (as in the case when T is even). Hence, we conclude that when T is even, the equation r_T+1 -(a+1)r_T +a = 0, r ≠ 1, -1 < a < 0 does not have any real root since f(r) tends to -∞ as r tends to -∞; when T is odd, the above equation has just one real negative root since f(r) tends to +∞ as r tends to -∞.

For the cases T = 1 to 4 it is possible to obtain explicit solutions in simple analytical form for the values of the eigenvalues r; for T > 4 the eigenvalues can be determined using numerical methods. Analytical procedures applied in the case T < 5 and extensive numerical analysis of the equation for any values of T > 4 have shown that (a) when T is odd, there is always one real eigenvalue, with |ri| < 1 and -1 < ri < 0, and (T-1)/2 pairs of complex conjugates r_j,k= A ± iB, with (A²+B²)^1/2< 1 for any initial gene frequency p₀; (b) when T is even, there are always only T/2 pairs of complex conjugates with the same properties as stated before. We therefore conclude that the set of equilibrium points p_e = f(p₀,T) is a stable one, with convergence always occurring to it for any initial value p₀in the open interval (0, 1) and for any integer value of T.

Description of the model in which there are fluctuations in the total number of individuals

In order to study the effect of an increase or decrease in population size, it is enough to introduce the additional parameter c: N₁p₁ =cN₀p₀into equation N₁p₁ =N₀p₀of the previous section. Consideration of this parameter leads straightforwardly the equation p_K+1 = 1/[2-p_K+(1-p_K-T)²/ (c_Tp_K-T)], K > T. This is the same equation obtained in the analysis of the constant population size model if c=1.

At equilibrium, p_K+1 = p_K = p_K-T = p_eand the above equation reduces to p_e(1-p_e)²(1-c_T)=0,c ≠ 1. When c ≠ 1, the only possible solutions for the foregoing equation are p_e= 0 and p_e =1, i.e., the normal allele is completely lost or fixed. Intuitively, it is easy to show that p_e = 0 when c < 1 (decreasing population size). Thus, for any K > T, the proportion of recessive genes introduced by larvae ending their diapause period will become increasingly greater as generations go by, in relation to the normal allele present in the adult population, the size of which is decreasing. The inverse, i.e., when c>1 (population increasing in size), is also obviously true since proportionally fewer recessive genes will be introduced into the adult population by emerging larvae. Clearly, then, p_etends to 1 as generations go by.

Brief remarks on a model with continuous generations

Let P(t), R(t), Q(t) and L(t) be the normalized proportions {P(t) + R(t) + Q(t) + L(t) = 1} of adults AA, Aa and aa, and larvae aa in generation t. Among adults, the normalized genotypic proportions are given by:

f(AA) = P' = P/(1-L)

f(Aa) = R' = R/(1-L)

f(aa) = Q' = Q/(1-L)

and the corresponding allelic frequencies f(A) and f(a) by

f(A) = p(t) = P' + R'/2

f(a) = q(t) = Q' + R'/2.

In the total population (adults + larvae), the allelic frequencies are given by

p(t)=P+R/2

and

q(t)=Q+L+ R/2, p(t) + q(t)= 1.

After a time interval Dt, the frequency of allele A among adults of the population is given by

p(t+Dt) = [N(t).p²(t) + N(t).2p(t)q(t)/2]/N(t+Dt) =

= {N(t).p(t)[p(t)+q(t)]}/N(t+Dt) =

= N(t).p(t)/N(t+Dt),

with the restriction N(t+Dt).p(t+Dt) = N(t).p(t) = N₀p₀= constant, as in the discrete model.

The expression for p(t+Dt) can be used for quantitative analysis (to be presented elsewhere) in the case of continuous generations. One expects that equilibrium values

for any fraction 0 < ε < 1 in T < T+ε < T+1 would interpolate between the equilibrium curves that represent the set of equilibrium points for the integer cases of T and T+1 generations in delayed larval development.

Acknowledgments

The subject of this paper was suggested to Paulo A. Otto in the 1980s by the late Prof. Luiz Edmundo de Magalhães, at the time a professor in the Department of Biology at USP. Subsequently, a congress note on this subject (Otto et al., 1990) was published with the co-authorship of Prof. Magalhães and a PhD student, Luiz Antonio Benedetti, now a professor in the Mathematics Department of the Federal University of Uberlândia in Minas Gerais. This work was partly funded by FAPESP.

Received: March 11, 2013; Accepted: July 7, 2013.

Associate Editor: Louis Bernard Klaczko

License information: This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

• Adell JC, Moya A and Botella LM (1988) Larval arrest in the development of the lygaeid bug Spilostethus pandurus Evol Biol 2:251-260.
• Bakker K (1959) Feeding period, growth and pupation in larvae of Drosophila melanogaster Entomol Exp Appl 2:171-186.
• Bakker K (1969) Selection for rate of growth and its influence on competitive ability of larvae of Drosophila melanogaster Netherl J Zool 19:541-595.
• Barker JSF and Podger RN (1970) Interspecific competition between Drosophila melanogaster and Drosophila simulans: Effects of larval density on viability, developmental period and adult body size. Ecology 51:170-189.
• Bialecki M, Shilton A, Fichtenberg C, Seagraves WA and Thummel CS (2002) Loss of the ecdysteroid-inducible E75A orphan nuclear receptor uncouples molting from metamorphosis in Drosophila Dev Cell 3:209-220.
• De Witt CT (1960) On competition. Versl Landbouwk Onderz 66:8-82.
• Dobzhansky T, Holz AM and Spassky B (1942) Genetics of natural populations. VIII. Concealed variability in the second and the fourth chromosomes of Drosophila pseudoobscura and its bearing on the problem of heterosis. Genetics 27:463-490.
• Dolezal T, Dolezelova E, Zurovec M and Bryant PJ (2005) A role for adenosine deaminase in Drosophila larval development. PloS Biol 3:e201.
• Ghosh A, McBrayer Z and O'Connor MB (2010) The Drosophila gap gene giant regulates ecdysone production through specification of the PTTH-producing neurons. Dev Biol 347:271-278.
• Huang SL, Singh M and Kojima KI (1971) A study of frequency-dependent selection observed in the esterase-6 locus of Drosophila melanogaster using a conditioned media method. Genetics 68:97-104.
• Mather K and Caligari PDS (1981) Competitive interactions in Drosophila melanogaster II. Measurement of competition. Heredity 46:239-244.
• McBrayer Z, Ono H, Shimell MJ, Parvy JP, Beckstead RB, Warren JT, Thummel CS, Dauphin-Villemant C, Gilbert LI and O'Connor MB (2007) Prothoracicotrophic hormone regulates developmental timing and body size in Drosophila Dev Cell 13:857-871.
• Meiklejohn CD, Holmbeck MA, Siddiq MA, Abt DN, Rand DM and Montooth KL (2013) An incompatibility between a mitochondrial tRNA and its nuclear-encoded tRNA synthetase compromises development and fitness in Drosophila PloS Genet 9:e1003238.
• Ménsua JL and Moya A (1983) Stopped development in overcrowded cultures of Drosophila melanogaster Heredity 51:347-352.
• Oldham S, Böhni R, Stocker H, Brogiolo W and Hafen E (2000) Genetic control of size in Drosophila Proc R Soc Lond B 355:945-952.
• Oliveira AK and Cordeiro AR (1981) Genetic control of developmental rate and a mother factor in D. melanogaster Rev Bras Biol 41:635-644.
• Oliveira AK, Monjelo LAS and Cordeiro AR (1991) A second chromosome major gene for developmental rate in Drosophila melanogaster Rev Bras Genet 14:953-965.
• Otto PA, Benedetti LA and Magalhães LE (1990) Efeitos do atraso de desenvolvimento larval sobre as frequências alélicas em uma população. Resumos do XXXVI Congresso Nacional de Genética, Caxambu, pp 167-167.
• Quinn L, Lin J, Cranna N, Lee JEA, Mitchell N and Hannan R (2012) Steroid hormones in Drosophila: How ecdysone coordinates developmental signalling with cell growth and division. In: Abduljabbar H (ed) Steroid Basic Science. Intech, Rijeka, pp 141-168.
• Robertson FV (1963) The ecological genetics of growth in Drosophila 6. The genetic correlation between the duration of the larval period and body size in relation to larval diet. Genet Res 4:74-92.
• Robertson FV (1964) The ecological genetics of growth in Drosophila 7. The role of canalization in the stability of growth relations. Genet Res 5:107-126.
• Sliter TJ and Gilbert LI (1992) Developmental arrest and ecdysteroid deficiency resulting from mutations at the dre4 locus of Drosophila Genetics 130:555-568.
• Tosic M and Ayala FJ (1981) Density-and frequency-dependent selection at the Mdh-2 locus in Drosophila pseudoobscura Genetics 97:679-701.

Publication Dates

• Publication in this collection
  23 Aug 2013
• Date of issue
  2013

History

• Received
  11 Mar 2013
• Accepted
  07 July 2013