Explicit solution of a Lotka-Sharpe-McKendrick system involving neutral delay differential equations using the r-Lambert W function

1.   Introduction

Delay differential equation (DDE) is a type of differential equation (DE) where the time derivatives at the present time are dependent on the solution and its derivatives at a previous time. A $k$-th order DDE takes the form

where $d_j = d_j(t, y(t))$ is called the delay satisfying $d_j \leq t$ for all times $t$ on the interval $[t_0, t_1]$, $j = 0, \ldots, k.$ A DDE is said to have discrete delays if the delays are constant. DDEs can have a single delay or multiple delays. A DDE is of retarded type (RDDE) when there is no time-delay in the derivative terms of the DDE. Its simplest form is

A DDE is of neutral type (NDDE) if there are delays in the derivatives or in the solution of the DDE. Its first-order form can be written as

Recall that the unique solution to an ordinary differential equation (ODE) is dependent on an initial condition. Analogously, the unique solution to a DDE is dependent on some initial function $\phi(t)$ defined on a previous time, say $[-h, 0]$, for some delay $h\in\mathbb{R^+}$. The function $\phi(t)$ is called the preshape or history function.

Asl and Ulsoy showed an approach for obtaining analytical solutions to systems of DDEs, particularly, first-order, linear, constant-coefficient RDDEs with constant delays using the Lambert $W$ function ^[1]. This was then extended by Yi to solve nonhomogeneous systems ^[2,3,4]. These results on DDEs are notable since there are few studies on explicit solutions to these types of differential equations.

The Lambert $W$ function, also known as the omega function, is defined to be the function $W(a)$ satisfying

If $a$ is a nonzero real number, then there are at most two real values of $W(a)$ satisfying (1.2). If $a$ is a complex number, then there is an infinite number of complex values of $W(a)$ satisfying (1.2). An in-depth discussion on the Lambert $W$ function, as well as the branch cut, asymptotes, and other properties, are found in ^[5].

A generalization of the Lambert $W$ function, introduced by Mező ^[6,7,8], is defined as a multi-valued inverse of the transcendental function

at $a$, where $a\in\mathbb{C}, \ a\neq0$ and $n, m\in\mathbb{N}$. We consider a simple form of the generalization,

with $n$ and $m$ equal to 1. Equation (1.3) can be transformed into the form as $xe^x+rx = a$, where $r\in\mathbb{R}$, and hence, is coined as the $r$-Lambert $W$ function. Therefore, the $r$-Lambert $W$ function is defined to be the function $W_r(a)$ satisfying the equation

In (1.4), we note the presence of the term $rW_r(a)$, which does not appear in the definition of the Lambert $W$ function in (1.2). If $a$ is a nonzero real number, then depending on the value of $r$, there can be one, two, or three real solutions of (1.4) with respect to $W_r(a)$. This means that there are at most three branches of the $r$-Lambert $W$ function, denoted $W_{(r, 0)}(x), W_{(r, -1)}(x)$, and $W_{(r, -2)}(x)$, having nonempty intersection with the real line. Extending to the complex plane, if $a$ is a complex number, then (1.4) has infinitely many complex solutions $W_r(a)$. This is similar to the solutions of (1.2) in the complex plane. However, Mező suggests that there are different structures of the branches of the $r$-Lambert $W$ function depending on the value of $r$ ^[8]. For a discussion on the solutions of (1.4), the branch structures, and other properties of the $r$-Lambert $W$ function, we refer the readers to ^[7,8]. The generalized Lambert $W$ function appears in the study of the solution to the Bose-Fermi mixtures ^[6] and on the asymptotic estimation of the Bell and generalized Bell numbers ^[9]. In ^[10], it was shown how the generalized Lambert $W$ function can be used to obtain a closed-form of the inverse of Langevin and Brillouin functions. These functions are used in several fields of physics. In ^[11], it was illustrated how the $r$-Lambert $W$ function can be used to characterize the equilibrium strategy in the symmetric two-player Hirshleifer contest. In ^[12], the $r$-Lambert function was used to show the asymptotic equivalence of the proximity operator of the logistic function and a composition of polynomial and exponential functions. For other applications of the generalized Lambert $W$ function in engineering, physics, and biology, we refer the readers to ^{[13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]}.

In ^[29], Jamilla et. al showed an approach for obtaining an analytical solution to first-order, linear, constant-coefficient NDDEs with constant delays using the $r$-Lambert $W$ function. The solution takes the form of a series whose terms are dependent on the parameters and the preshape function of the NDDE. In this paper, we adopt an age-structured population model described by a system of ODEs and NDDEs. We utilize the $r$-Lambert $W$ function to provide an explicit solution to the model using the approach proposed in ^[29]. One advantage of this method is that its form is similar to the form of the solution of ODEs ^[1]. Moreover, as stated above, the terms in the obtained solution are dependent on the parameters of the NDDE, and hence, one can analyze the effect on the solution of varying parameter values.

The outline of this paper is as follows. We discuss in section 2 the Lotka-Sharpe-McKendrick age-structured population model presented in ^[30]. We show how the model is reduced to a system of ODEs and NDDEs, along with the necessary conditions and relationships of the parameters of the reduced and original model. In section 3, we solve the reduced system using the $r$-Lambert $W$ function. An approach to numerically approximate the solution obtained in section 3 is detailed in section 4. Because the true solution of the model is not available, we compare our method to the MATLAB built-in solvers. We test our method by varying essential parameters and initial conditions, which we illustrate in three examples. The summary, conclusion, and future works are in section 5.

2.   Age-structured model

Bocharov and Hadeler in ^[30] considered a linear Lotka-Sharpe-McKendrick system of an age-structured population given by

where $n(t, a)$ is a non-normalized age distribution of the population at time $t$, $a$ is the chronological age, $\mu(a)$ is the mortality rate, and $b(a)$ is the fertility rate. Both $\mu(a)$ and $b(a)$ are functions dependent on age, and are assumed to be sums of step functions and delta peaks of the following forms:

where

$\delta_\tau(a)$ is the Dirac delta function, and $\mu_k$ and $b_k$, $k = 0, 1, 2$, represent death and birth rates, respectively. We assume that $\mu_k$ and $b_k$ are nonnegative. The coefficients $\mu_k$ and $b_k$ are chosen so that they mimic the qualitative behavior of real populations such as non-constant mortality, reproductive window, etc. ^[30]. In Figure 1, we see three coefficients $b_0, \ b_1, \ b_2$ that constitute the fertility rate $b(a)$. When $a < \tau$, then $b(a) = b_0$, while if $a > \tau$, $b(a) = b_1$. Notice that $b(a)$ has a peak at $a = \tau$. Similarly, the coefficients $\mu_0$ and $\mu_1$ appear in $\mu(a)$, and the graph of $\mu(a)$ has a jump at $a = \tau$. The functions for the mortality and fertility rates are chosen such that the rates are different when $a\leq \tau$ and when $a > \tau$. Based on this, the population is divided into two classes, the juveniles ($a\leq \tau$) and the adults ($a > \tau$), each with corresponding mortality and fertility rates.

The following new variables are then introduced:

where $U(t)$ is the population of the juveniles and $V(t)$ is the population of the adults at time $t$. We assume that the juveniles do not reproduce ($b_0 = 0$). For example, the nymph or larval stage of most insects lacks functional reproductive organs ^[31]. We also assume that there is no sudden increase in the death rate at the age $\tau$ ($\mu_2 = 0$). Using the functions $U(t)$ and $V(t)$, the partial differential Eq (2.1) is transformed into a system of the form

with initial conditions

dependent on $n_0(a)$.

The parameters $a_i$ and $c_j$ are related to the parameters $\mu_k$ and $b_k$ in the following manner:

where $i = 1, 2, 3$ and $j = 0, 1, 2, 3$. Moreover, from ^[30], it is necessary that $a_i\geq0$ for all $i$ and that $a_1\geq a_2a_3$.

We see that the partial differential equation in (2.1) is reduced to a system involving ODEs and NDDEs. Throughout this paper, we consider the age-structured population model (2.6)–(2.9). We present an explicit solution to this model and illustrate numerically its solution using the $r$-Lambert $W$ function.

3.   Solution to the model

In this section, we consider system (2.6)–(2.7), along with the initial conditions (2.8) and (2.9). First, we determine the solutions $U(t)$ and $V(t)$ to the model for $t\in\left[0, \tau\right]$. The solutions on $\left[0, \tau\right]$ are the preshape functions necessary to find the solutions $U(t)$ and $V(t)$ for $t > \tau$.

3.1.   Solutions for $t\in\left[0, \tau\right]$

We begin by solving for $U(t)$ and $V(t)$, $t\in[0, \tau]$, from the ODEs in (2.6) and (2.7). Note that $V(t)$ is independent of $U(t)$, and thus, can be solved separately. Applying techniques in solving ODEs, we get

where $C_1$ is an arbitrary constant. Next, we solve for $C_1$ using the initial condition $V(0)$. We let

Substituting (3.2) to (3.1) and solving for $C_1$ at $t = 0$, we get

Therefore,

Next, we solve for $U(t)$ for $t\in\left[0, \tau\right]$. We plug-in the expression for $V(t)$ from (3.3) to the first equation in (2.6). The resulting equation is a non-homogeneous ODE, and its solution is given by

where $C_2$ is an arbitrary constant. We then solve for $C_2$ by letting

and

Using the initial condition $U(0)$, we get

Thus,

3.2.   Solutions for $t > \tau$

We now solve system (2.6)–(2.9) for $t > \tau$. Since $V(t)$ is independent of $U(t)$, we solve $V(t)$ first following the methodology in ^[29]. Suppose that the solution $V(t)$ to the second equation in (2.7) is of the form $e^{\lambda t}$. Plugging in $V(t) = e^{\lambda t}$ into the NDDE, we get

Multiplying both sides of the equation by $\tau e^{\tau(\lambda + a_3)-\lambda t}$ yields

Now, adding a constant term $-a_2\tau a_3e^{\tau a_3}$ to both sides of the equation and rearranging terms, we obtain

which is of the form $xe^x + rx = a$, with $x = \tau(\lambda + a_3),$ $a = a_1\tau e^{\tau a_3}-a_2\tau a_3e^{\tau a_3}$, and $r = -a_2e^{\tau a_3}$. Using the $r$-Lambert $W$ function, we can write $\lambda$ as

where $W_r$ is defined to be solution of the $r$-Lambert $W$ function. Noting that there is an infinite number of solutions to the $r$-Lambert $W$ function, we express $V(t)$ as

where

$W_{(r, k)}$ is defined to be the $kth$ branch of the $r$-Lambert $W$ function with $r = -a_2e^{\tau a_3}$, and $C_k$ is dependent on $V(t)$ for $t\in[0, \tau]$.

Next, we solve for $U(t)$ for $t > \tau.$ From (3.7), we obtain

and

Substituting these expressions to (2.6), we get

Suppose that a solution to the associated homogeneous DE $\dot{U}(t) = c_0 U(t)$ is of the form $U(t) = e^{\lambda t}$. Then $\dot{U}(t) = \lambda e^{\lambda t}$, and hence, $\lambda = c_0$. Thus, the complementary function for (3.9) is

for some constant $A$. To find a particular solution $U_p(t)$ to (3.9), suppose that

Then

From (3.9), we have

Rearranging the terms, we get

where

Comparing coefficients on both sides of (3.11), it follows that for each $k$,

and

Then

and

Now, we solve for $A$ in (3.10). Since we are solving for $U(t)$ for $t > \tau$, then it follows from (3.12) that the initial condition is given by

Hence,

Finally, the solution to system (2.6)–(2.9) is given by

where

and $C_k$ is dependent on $V(t)$ on $[0, \tau]$.

4.   Numerical illustration

Note that the explicit solution computed in the previous section for $t > \tau$ involves an infinite series, and thus, cannot be exactly solved. We can only approximate the solution by limiting the summations, i.e., limiting the number of solutions to the $r$-Lambert $W$ function being considered. For $N\in\mathbb{N}$ sufficiently large, we set

on the interval $[\tau, T]$, where $T > \tau$. Also, note that the solution $V(t)$ on $(\tau, \infty)$ is dependent on a history function $\phi(t)$, which is $V(t)$ on the interval $[0, \tau]$. To solve for $V(t)$ numerically on MATLAB, we need an approximate for the coefficient vector $C_k$. For this, we use the procedure discussed in ^[29]. The history function $\phi(t)$ can be approximated as follows:

where $\zeta_k(t) = e^{\lambda_kt}$, with the same $\lambda_k$ as in (3.8). Partitioning the interval $[0, \tau]$ to get $2N+1$ points, equation (4.2) becomes

Solving the system above gives the coefficients $C_k, \, k\in\{-N, \ldots, N\}$ in (4.1).

We illustrate in the following three examples the solution obtained in section 3 using the approximation (4.1). In the first two examples, we vary the initial age distribution $n_0 (a)$. In the third example, we investigate under certain conditions the long-term behavior of the solution. The parameter $t$ represents time in months. Because we do not know the true solution to system (2.6)–(2.9), we compare our obtained solution with the solution using the MATLAB built-in functions $\texttt{ode23}$ and $\texttt{ddensd}$.

In the next two examples, we use the parameter values shown in Table 1.

We set $T = 144$ and compute for the following parameters:

Example 1. We set the initial age distribution to be uniform, that is,

where $a_{max}$ is the maximum age in the population. Thus, the initial conditions are

For $t\in[0, 36],$ the solutions $V(t)$ and $U(t)$ are given by

with

Moreover, the approximate solutions $U(t)$ and $V(t)$ on the interval [36, 144] are given by

and

where

and $C_k$ is a column vector dependent on $\phi(t) = V(t)$ for $t\in[0, 36]$.

We evaluate the explicit solutions $U(t)$ and $V(t)$ on the interval $[0, 144].$ Note that we approximate the solutions $U(t)$ and $V(t)$ on the interval [36, 144] using $N = 50$. We also compute for the solution using the built-in functions $\texttt{ode23}$ and $\texttt{ddensd}$ from MATLAB with stepsize $s = 0.01$. We then plot our proposed solutions, $U(t)$ and $V(t)$, together with the solutions obtained using the MATLAB built-in functions, denoted by $U_M(t)$ and $V_M(t)$, for 30 equally spaced points on $[0, 144]$. Figures 2a and 2b show the plots of these results.

Notice that there is a decline in the juvenile population for the first 36 months, compared to the steady increase in the adult population. Moreover, there is an almost steady increase for both populations for $t > 36.$ Also, at around $t = 76$ months, the total population is doubled, and at $t = 144$ months, the total population is quadrupled.

Comparing our solutions $U(t)$ and $V(t)$, with the numerical solutions $U_M(t)$ and $V_M(t)$ using the built-in functions at 30 different points, we see that the two solvers obtain similar results. The relative differences in approximation are

Going back to system (2.6)–(2.9) and the parameter values (2.10), notice that the absence of the parameter $b_2$ will transform the NDDEs into RDDEs. Thus, we consider $b_2$, together with the delay $\tau$, as the two most significant parameters. To show that our method is robust, we implement our method by varying the parameter values of $b_2$ and $\tau$. We then compare our results with the MATLAB built-in functions using the following measure of relative difference

where $t_i, \in[0, 144], i = 1, ..., 30$. We use the same values in Table 1 for the parameters $b_1$, $\mu_0$, and $\mu_1$, and choose $b_2$ and $\tau$ such that $b_2\in[0.0013, 0.0045]$ and $\tau\in[24, 48]$. In particular, we evaluate $R(b_2, \tau)$ at 50 different values of $b_2$ and $\tau$, and create a heat map to compare the obtained results. Figures 3a and 3b show the values of $R$ in (4.4) for different values of the parameters $b_2$ and $\tau$ when $N = 10$ and $N = 50$, respectively.

Notice that the relative difference when $N = 50$ is smaller than when $N = 10$, for all $b_2$ and $\tau$. Also, the minimum value of $R$ is $4.75501\times10^{-8}$ when $N = 10$, while the maximum value is $1.83404\times10^{-5}.$ When $N = 50$, the maximum value of $R$ is $6.96986\times10^{-8}$, and the minimum value is $3.75354\times10^{-10}$. This shows that when $b_2$ varies between $[0.0013, 0.0045]$ and $\tau$ varies between $[24, 48]$, our proposed method can give a good approximate solution.

Example 2. We also consider a normal distribution as the initial age distribution, given by

where $\mu$ and $\sigma$ are the mean and standard deviation, respectively. The values of $\mu$ and $\sigma$ are set to 25 and 5, respectively. This distribution follows a bell-shape structure and indicates that the age $a$ with the highest number of initial population is 25. The plot of $n_0(a)$ is illustrated in Figure 4.

Again, we set the maximum age to be 50. Solving for the initial conditions, we have

We let $t = \dfrac{a-25}{5\sqrt{2}}$. Then we have $dt = \dfrac{da}{5\sqrt{2}}.$ Transforming (4.8) and (4.9), we have

and

Note that the addends above for $U(0)$ and $V(0)$ are similar to the form of the error function defined as

The error function arises when integrating the normal distribution. Substituting $\text{erf} (x)$ to $U(0)$ and $V(0)$, we get

and

For $t\in[0, 36],$ the solutions $V(t)$ and $U(t)$ are given by

with

and

Approximate solutions $U(t)$ and $V(t)$ on the interval $[36, 144]$ are given by

and

where

and $C_k$ is a column vector dependent on $\phi(t) = V(t)$ for $t\in[0, 36]$.

Again, we evaluate the solutions $U(t)$ and $V(t)$ on the interval $[0, 144]$ using the parameter values in Table 1. We approximate the solutions on the interval $[36, 144]$ with $N = 50$. We also compute for the solutions using the MATLAB built-in functions $\texttt{ode23}$ and $\texttt{ddensd}$ with stepsize $s = 0.01$. We then compare our solutions $U(t)$ and $V(t)$, with the numerical solutions $U_M(t)$ and $V_M(t)$ obtained using the built-in functions, for 30 equally spaced points in $[0, 144]$. We see in Figures 5a and 5b that the graphs of $U(t)$ and $V(t)$ are close to the graphs of $U_M(t)$ and $V_M(t)$, respectively. The relative differences in approximation are

Similar to the previous example, there is a decline in the population of the juveniles during the first 36 months, while there is an increase in the population of the adults. We can also see that for $t > 36,$ both populations increase. In fact, at $t = 95$, the total population is twice the initial total population. It takes more time to double the population in this example, compared to Example 1.

We also compare the difference in approximation of our proposed method and the MATLAB built-in functions. We calculate $R(b_2, \tau)$ in (4.7) at 50 different values of $b_2$ and $\tau$, and plot the results when $N = 10$ and $N = 50.$ Figures 6a and 6b show the results. As expected, the values are smaller when $N = 50$ than when $N = 10$. When $N = 10$, the maximum value of $R$ is $4.51858\times10^{-7}$, while the minimum value is $1.79839\times10^{-9}.$ When $N = 50$, the maximum value of $R$ is $6.96986\times10^{-8}$ and the minimum value is $3.75354\times10^{-10}.$ The results show that our proposed explicit formulation can be used to numerically solve (2.6)–(2.9).

Example 3. In this example, we consider the model studied by Gourley and Kuang in ^[34]. The model is a system of ODEs and NDDEs given by

where $u_i$ and $u_m$ represent the immature (juvenile) and mature (adult) members of the population, respectively, and $u_0(a)$ represents the initial age distribution of the corresponding PDE system, similar to system (2.1). Moreover, they considered a similar birth rate function (2.5) and set $\mu$ as the constant linear death rate for the juveniles and $d(u_m(t))$ as the adult mortality function satisfying

System (4.8)–(4.9) can be used to describe the population of insects with a very long larval stage and a short adult stage allotted for mating ^[34]. One example of such insects is the 17-year periodical cicadas, particularly the species Magicicada septendecim, M. cassini, and M. septendecula, which spend years of their life as nymphs, feeding underground on plant root xylems and resurfacing as adults only to mate for a few weeks ^[35]. Another example is the marine midges of genera Pontomyia and Clunio, which spend a month as benthic larvae and a few hours as adults ^[36]. Mayflies of order Ephemeroptera are another example, since they generally live from 3–4 weeks to more than 2 years as nymphs, and from a few hours to a few days as adults ^[37]. Note that the juveniles of these insects are not capable of reproduction and hence, we set the juvenile reproduction to be 0 ($b_0 = 0$).

Assuming that the adult mortality function $d(u_m)$ is linear, i.e.,

for some $k > 0$, system (4.14)-(4.15) is reduced to

A theorem on the existence of a solution to system (4.14)-(4.15) is proven in ^[38], and this also holds for the reduced system. Note that system (4.17)–(4.18) is composed of linear ODEs and NDDEs and thus, using our method, a solution to this system is given by

where

and $C_j$ is dependent on $u_m(t)$ on $[0, \tau]$. Table 2 lists the parameter values used in the simulations. Furthermore, we set $k = 0.5$ and the initial age distribution to be

The chosen parameter values for $b_0,$ $b_1$ and $b_2$ in Table 2 imply that reproduction occurs at a very short period when the insect reaches adulthood. Figures 7a and 7b show the plots of the immature and mature populations (labeled $u_i$ and $u_m$, respectively) using (4.19)–(4.20) on the interval $[0, 30]$ with $N = 50$, and using the MATLAB built-in functions $\texttt{ode23}$ and $\texttt{ddensd}$ with $s = 0.01$ (labeled $u_i^*$ and $u_m^*$, respectively).

The relative differences in approximation of the solutions for 30 equally spaced points in $[0, 30]$ are $R_{u_i} = 3.7641 \times 10^{-07}$ and $R_{u_m} = 3.7869 \times 10^{-06}$, which imply that $u_i(t)$ and $u_m(t)$ are close to $u_i^*(t)$ and $u_m^*(t)$, respectively. Moreover, we observe in Figure 7b that as $t\to\infty$, $u_m(t)$ converges to 0. Also, since the juvenile population is dependent on the adult population, $u_i(t)$ also converges to 0 (Figure 7a). This convergence agrees with Theorem 2 of ^[34], which states that if $b_0 = b_1 = 0$ and $b_2e^{\mu t} < 1,$ and $d(u_m(t))$ is a continuous strictly monotonic increasing function of $u_m$ satisfying $d(0) = 0$, then $u_m\to0$ as $t\to\infty.$

Next, we consider the case when $b_1\ne0$. We verify the asymptotic behavior of the resulting system using Theorem 3 in ^[34], which states that if $b_0 = 0,$ $b_1 > 0,$

and $d(u_m)$ is a continuous strictly monotonic increasing function of $u_m$ satisfying $d(0) = 0,$ then $u_m\to0$ as $t\to\infty.$ Using our assumption on $d(u_m)$ and simplifying the inequality above, we infer that if we set

then $u_m\to0.$ Consequently, using the parameter values for $b_0, b_2, \mu$ and $\tau$ in Table 2 and setting $b_1 = 2$, then $u_m\to0$ if $k > 9.1296$. Figures 8a and 8b show our solutions $u_i(t)$ and $u_m(t)$ with $k = 13$ and $N = 200$, in comparison with the solutions $u_i^*(t)$ and $u_m^*(t)$ using the MATLAB built-in functions with $s = 0.01$. We see in Figures 8a and 8b that $u_i\to0$ and $u_m\to0$ as $t\to\infty.$ Moreover, the relative differences of the solutions for 30 equally spaced points on $[0, 50]$ are $R_{u_i} = 8.0737\times10^{-07}$ and $R_{u_m} = 1.1265\times10^{-05}$. This implies that our solution is close to the solution obtained using the MATLAB built-in solvers.

In this example, we have shown how our proposed numerical method can be used to investigate the asymptotic behavior of the solution.

5.   Conclusion and future works

We have presented an explicit solution to an age-structured model. By assuming that the juveniles do not reproduce and that there is no peak in the death rate when juveniles become adults, the age-structured model is reduced to a system of ODEs and NDDEs. We have shown that the arising NDDE can be solved using a generalization of the Lambert $W$ function. An explicit solution in the form of an infinite series was obtained. The terms of this series depend on the parameters of the NDDE. Hence, one can identify the effects of changing parameter values in the solution. We tested our approach numerically and compared our computed solution with that of the MATLAB built-in solvers.

For future work, one may study the solution of the problem when there is no assumption on the non-reproduction of juveniles and no peak in the death rate when juveniles become adults. In Example 3, we have shown that our proposed numerical approach can be used to study the long-term behavior of the solution. However, an in-depth theoretical analysis of the convergence of the solution of the model using our proposed explicit solution is not trivial. This would require a thorough analysis of the $r$-Lambert $W$ function and would demand a separate study. The explicit solution derived in this work relies heavily on the structure of the solutions of (1.4). To deal with a nonlinear model, one needs to study the structure of the solutions of (1.4) that is not only dependent on the constant $a$ but also on time $t$. This is an exciting direction for future research.

Acknowledgments

The authors acknowledge the Office of the Chancellor of the University of the Philippines, through the Office of Vice Chancellor for Research and Development, for funding support through the Outright Research Grant.

Conflict of interest

All authors declare no conflicts of interest in this paper.