Bifurcation analysis of a mathematical model of microalgae growth under the influence of sunlight BIOMATH https://biomath.math.bas.bg/biomath/index.php/biomath B f Biomath Forum ORIGINAL ARTICLE Bifurcation analysis of a mathematical model of microalgae growth under the influence of sunlight Lingga Sanjaya Putra Mahardhika, Fajar Adi-Kusumo∗, Dwi Ertiningsih Department of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Gadjah Mada, Yogyakarta, Indonesia linggasanjaya2018@mail.ugm.ac.id 0000-0001-5063-2848 f adikusumo@ugm.ac.id 0000-0002-1643-4466 dwi ertiningsih@ugm.ac.id 0000-0002-4529-8440 Received: July 9, 2022, Accepted: January 30, 2023, Published: April 28, 2023 Abstract: In this paper is considered a microalgae growth model under the influence of sunlight. The model is a two-dimensional system of the first order Ordinary Differential Equations (ODE) with a ten-dimensional pa- rameter space. For this model, we study the existence of equilibrium points and their stability, and determine a bifurcation of the system when the value of some parameters is varied. The Lambert ω function is used to calculate equilibrium points and apply the linearization technique to provide their stabilities. By varying the value of some parameters numerically, we found a transcritical bifurcation of the system and show stability regions of the equilibrium points in parameter diagrams. The bifurcation shows that the microalgae have a minimum sustainable nutrition supply and have a minimum light intensity that plays an important role for survival in a chemostat which has a certain depth. The results can be used to design a chemostat in optimizing the growth of microalgae. Keywords: Microalgae growth model, Quota cell, Pa- rameter diagram, Bifurcation I. INTRODUCTION At the moment, fossil fuels still generate about 80% of the demand of global energy. This demand increases along with the increase in population, because each individual needs a means of transportation to carry out activities and move to other places. However, the extensive use of fossil fuels plays an important role for global climate change, environmental pollution, and problems in health [1]. Most scientists are looking for new types of energy. The most interesting renewable energy that is expected to have an important role in the future global energy structure is biofuels. Biodiesel is one of the biofuels that is recognized as an ideal carrier of renewable energy which has the potential to become a primary energy source. There are several candidates of plants for biodiesel production, but most of them grow slowly, contain little vegetable oil for biodiesel, and require large areas of land to be grown. Therefore, they are considered inefficient for biodiesel production. Microalgae is a microorganism that has an ability to convert solar en- ergy to chemical energy through the fixation of carbon dioxide (CO2). The growth of microalgae is relatively fast, so it can be considered as a valuable source for making biodiesel [1]. Microalgae can form their energy by photosynthesis to support their life needs. It carries out photosynthesis with the help of sunlight and absorbs carbon dioxide and nutrients around it to form the various substances it needs. The result of photosynthesis is glucose that can be stored in cells as vegetable fats. Copyright: © 2023 Lingga Sanjaya Putra Mahardhika, Fajar Adi-Kusumo, Dwi Ertiningsih. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. *Corresponding author Citation: Lingga Sanjaya Putra Mahardhika, Fajar Adi-Kusumo, Dwi Ertiningsih, Bifurcation analysis of a mathematical model of microalgae growth under the influence of sunlight, Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 1/9 https://biomath.math.bas.bg/biomath/index.php/biomath mailto:linggasanjaya2018@mail.ugm.ac.id https://orcid.org/0000-0001-5063-2848 mailto:f_adikusumo@ugm.ac.id https://orcid.org/0000-0002-1643-4466 mailto:dwi_ertiningsih@ugm.ac.id https://orcid.org/0000-0002-4529-8440 https://creativecommons.org/licenses/by/4.0/ https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight Biofuel is mostly made from vegetable oil or animal fats [2]. Therefore, microalgae can play an important role to reduce carbon dioxide levels in the atmosphere and become a candidate for renewable energy sources to replace fossil fuels [3]. Furthermore, the potential of biofuel products from the microalgae is expected to be abundant [4]. The theoretical investigation of microbial growth and nutrient dynamics in chemostats was pioneered by J. Monod [5]. The model proposed that extracellular nu- trient concentration is a determinant of cellular growth. However, microorganisms have nutrient reserves in cells that are used for growth and development, which causes the growth rate of microorganisms depend on the nutrient reserves in cells. Thus the Monod model is a weak growth predictor. Furthermore, M.R. Droop [6], purposed a new empirical model as an update of the Monod model which is known as the cell quota model with a specific growth rate. Droop’s model is a function of cell quota (intracellular density) of limiting nutrients. The Monod and Droop models are the basis of the general microorganism growth model. In fact, mi- croalgae not only absorb nutrients but also perform photosynthesis. The growth model of microalgae is mostly based on the Droop model and then modified with the addition of the sunlight factor. However, the model is limited by the conditions of low nutrient levels. Since light intensity which transmits a medium can be explained with Lambert-Beer law [7], the maximum growth of microalgae under influence of light can be formulated by Lambert-Beer law. Furthermore, the microalgae growth model – this is a new derivation after the Lambert-Beer law is considered – contains the exponential term. There are several problems in the mathematical models of microalgae growth. One of which is the model with the influence of sunlight that contains an exponential form, so it is quite difficult to carry out stability analysis and bifurcation analysis on the model. Our model is motivated by the one in [5], where the mortality rate of the system in [5] represents chemostat dilution. However, based on the fact that microalgae have a natural mortality rate at room temperature, see [8], we add the microalgae natural mortality rate in our system. Since this model contains an exponential term, we apply the Lambert ω function to determine equilibrium points (see [9]). Furthermore, in the fourth section, we analyze the stability of the equilibrium point by linearization technique (see [10, 11]) to provide the stability of the equilibrium points. Since the stability conditions depend on the parameters, then we can di- vide the parameter space into several regions. Lastly, the occurrence of bifurcation is investigated by analyzing those regions numerically (see the bifurcation theory in [12–14]) and are established bifurcation parameters with respect to four varying parameters. II. MICROALGAE GROWTH MODEL A. Microalgae Microalgae are autotrophs that produce their food through the process of photosynthesis [15]. Photosyn- thesis is the process of converting inorganic compounds CO2 and H2O into glucose under sunlight. Solar energy is used in this process [16, 17]. Nutrients in cell quotas depend on nutrients in the environment, while microalgae growth depends on pho- tosynthesis and nutrient levels in cell quotas. Microal- gae require large amounts of macro and micro nutrients for their growth. There is a linear relationship between nutrient density and biomass, because microalgae ab- sorb nutrients from the environment and then store them in the cell quotas [1]. Microalgae assimilate various organic carbon and inorganic sources as nutrients (such as glucose, acetate, nitrogen, and phosphor) for their growth [18]. Light is used by microalgae to break down CO2 into glucose in the process of photosynthesis. To make an efficient microalgae culture for microalgae growth, it is necessary to optimize the intensity of light that enters the culture. In indoor and outdoor microalgae cultivation systems, light source and light intensity are important factors that affect the phototrophic growth performance of microalgae [1]. Light can be transmitted by a medium that absorbs the light intensity. The light intensity that is absorbed depends on the concentration of liquid and the thickness of the medium [19]. Microalgae have a level of turbidity in the water so that they can block the light intensity. Thus, if the microalgae content in the chemostat is denser, less light can enter. So, less light can be used by microalgae. Several models in the literature have been developed to explain microalgae growth under influence of light, especially in monocultures [20]. B. Microalgae growth models under influence of sun- light We propose a new microalgae growth model under the sunlight influence which is motivated by the one in [5]. In our model, we consider the microalgae natural mortality rate at room temperature m, not due to harvesting D. This could explain more the natural Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 2/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight behavior of microalgae at room temperature. The model is as follows: dA(t) dt = A(t) ( lN(t)βkI0e −(σ+kA(t))z lN(t) + qβkI0e−(σ+kA(t))z −m ) (1) dN(t) dt = D ( Nc −N(t) ) − lN(t)A(t), (2) where A(t),N(t),σ,k,I0,Nc and z represent the mi- croalgae content, nutrient content, turbidity of the wa- ter, the coefficient of turbidity due to the presence of microalgae concentration, light intensity, stock of nutrient, and depth of the medium respectively. Note: Nc, A(t), q, and N(t) were normalized by the density of water (1 g cm−3) to obtain relative density (n.d. – non-dimensional) in the model (see [5]). The chemostat contains N(t) (in g cm−3) extra- cellular nutrients for microalgae. Microalgae A(t) (in g cm−3) in the chemostat is always given new nutrients Nc in the form of dilution D to replace the old nutrients in the chemostat. Furthermore, the new nutrients enter the chemostat with the rate DNc and the old nutrients are pushed out through the pipe with the rate DN(t). According to Lambert-Beer law, the light intensity can be absorbed by microalgae at the rate kA(t)z. The thickness of the water causes the light intensity to decrease with a decreasing rate σz. By using Lambert- Beer law equation, the light intensity can be formulated as I = I0e−(σ+kA(t))z. Furthermore, microalgae need sunlight to carry out photosynthesis and then the prod- ucts of photosynthesis are used to grow. We assume the maximum microalgae growth rate has a linear relationship with the light intensity. The relationship between the maximum microalgae growth rate and the light intensity can be formulated as µmax = βkI, where β represents a constant the ratio between the maximum microalgae growth rate and the light intensity. Microalgae absorb nutrients from outside the cells into the quota cell at the rate lN(t). According to the Droop quota cell model, the nutrients in the quota cells help the microalgae to grow. Then the nutrient in the quota cell can be formulated Q = lN(t) µmax + q. By using the Droop quota cell model, the microalgae growth rate can be formulated as µ = µmax(1 − qQ), where q and µ represent the minimum quota cell for microalgae to grow and the growth rate of microalgae, respectively. Furthermore, the microalgae content is increased by µA(t). Theorem 1 guarantees a non-negative model solution. Next, Theorem 2 guarantees that the model solution is bounded. Theorem 1. Given A(0) and N(0) are non-negative. If all parameters in the model are positive, then A(t) and N(t) are non-negative. Proof: Will be proved by contradiction. Assume that the solution A(t) and N(t) always decrease and then will be negative. Let say that as follows. For 0 < t′ < t0, then A(t′) > 0 and N(t′) > 0. At the time t0, the variables are A(t0) = 0 and N(t0) = 0, and its derivatives are A′(t0) < 0 and N′(t0) < 0. Substituting t0 into Equation (1) and (2), then we get the following: dA(t0) dt = A(t0) ( lN(t0)βkI0e −(σ+kA(t0))z lN(t0) + qβkI0e−(σ+kA(t0))z −m ) = 0, dN(t0) dt = D ( Nc −N(t0) ) − lN(t0)A(t0) = DNc > 0. So that, there is a contradiction. Hence, the solution A(t) and N(t) will not be negative. Theorem 2. Given A(0) and N(0) are non-negative. If all parameters in the model are positive, then the solution A(t) and N(t) are bounded. Proof: Will be proved by contradiction. Assume that the solution A(t) and N(t) are unbounded. Hence, A(t) and N(t) always increase. For N(t) will be: N(t) > Nc. Let say that as follows. For 0 < t′ < t0, at the time t′: N(t′) < Nc. At the time t0, the variables are: N(t0) = Nc. Because the solution N(t) always increase, then dN(t0) dt > 0. Substituting t0 into Equation (2), then we get the following: dN(t0) dt = D (Nc −N(t0)) − lN(t0)A(t0) = −lNcA(t0) < 0. So that, there is a contradiction. Hence, the solution N(t) is bounded. Furthermore, N(t) will not be greater than Nc. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 3/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight Define r1 = lβkI0, r2 = qβkI0. Now for the solution A(t), considered below: dA(t) dt = A(t) ( lN(t)βkI0e −(σ+kA(t))z lN(t) + qβkI0e−(σ+kA(t))z −m ) = A(t) ( N(t)r1 lN(t)e(σ+kA(t))z + r2 −m ) = A(t) ( N(t) ( r1 −mle(σ+kA(t))z ) −mr2 h1 ) , for h1 = lN(t)e(σ+kA(t))z + r2. Since A(t) always increase, then dA(t)/dt > 0 for all t > 0. It means, N(t) ( r1 −mle(σ+kA(t))z ) −mr2 h1 > 0, for all t > 0. Since h1 > 0, then: N(t) ( r1 −mle(σ+kA(t))z ) > mr2. Since N(t) and mr2 are non-negative, then we get as follows: r1 −mle(σ+kA(t))z > 0. Furthermore, note that A(t) always increase. Let say that as follows. For 0 < ts < t̄, such that at the time ts we have: r1 −mle(σ+kA(ts))z > 0. Because the solution A(t) always increase, at the time t̄ we get: r1 −mle(σ+kA(t̄))z = 0, and dA(t̄)/dt > 0. Next step, substituting t̄ into Equation (1), then we get: dA(t̄) dt = A(t̄) N(t̄) ( r1 −mle(σ+kA(t̄))z ) −mr2 h1 = −A(t̄) mr2 h1 < 0. So that there is a contradiction. Hence, the solution A(t) is bounded. III. DETERMINE EQUILIBRIUM POINTS In this section we will discuss about how to de- termine equilibrium points of the model. The Lam- bert ω function is used to determine the equilibrium point. Equilibrium conditions are dA(t)/dt = 0 and dN(t)/dt = 0, so that the following equation is obtained: 0 = A(t) ( lN(t)βkI0e −(σ+kA(t))z lN(t) + qβkI0e−(σ+kA(t))z −m ) , (3) 0 = D (Nc −N(t)) − lN(t)A(t). (4) From the equation (3) we divide it into two cases, they are A = 0 or A 6= 0. For A = 0 obtained below: 0 = D (Nc −N(t)) − lN(t)A(t) ⇒ N(t) = Nc. So, the first equilibrium point is (A∗1,N ∗ 1 ) = (0,Nc). For A 6= 0 in the nutritional equation dN(t)/dt = 0, it is obtained as follows: 0 = D (Nc −N(t)) − lN(t)A(t) ⇒ N(t) = DNc D + lA(t) . (5) Furthermore, because A 6= 0, in order to satisfy the equilibrium condition, it is obtained as follows: N(t) ( r1 −mle(σ+kA(t))z ) −mr2 = 0. (6) Substituting Equation (5) into Equation (6) yields: DNc ( r1 −mle(σ+kA(t))z ) −mr2(D + lA(t)) D + lA(t) = 0. (7) Prior to further discussion, note the following coeffi- cients. Recall that r1 = lβkI0, r2 = qβkI0. Define a = mlDNce σz, b = lmr2, c = Dmr2 −DNcr1, and d = kz. So that to meet the equilibrium conditions, (7) yields: aedA(t) + bA(t) + c = 0. (8) Furthermore, Equation (8) is an exponential equation, using the Lambert ω function (for example see [21]), the root of the equation is as follows: A∗ = − bω ( ade−cd/b b ) + cd bd , (9) where ω ( ade−cd/b b ) is the Lambert ω function that is executed on ade −cd/b b . By substituting Equation (9) into Equation (5), we get the following: N∗ = DNcbd Dbd− l ( bω(ade −cd/b b ) + cd ). (10) So, the second equilibrium point (A∗2,N ∗ 2 ) is:( − bω(ade −cd/b b ) + cd bd , DNcbd Dbd− l ( bω(ade −cd/b b ) + cd ) ) . Furthermore, the existence of these two equilibrium points will be investigated. In this model, the equilib- rium point exist, if A∗,N∗ ≥ 0. Lemma 3. If all parameters in the model are positive, then the equilibrium point (A∗1,N ∗ 1 ) = (0,Nc) exists. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 4/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight Proof: It is known that all parameters are greater than 0, then Nc > 0. Since A∗1 = 0 and N ∗ 1 = Nc > 0, the equilibrium point (A∗1,N ∗ 1 ) = (0,Nc) exists. Lemma 4. If all parameters in the model are positive and a ≤ −c, then the equilibrium point (A∗2,N∗2 ) =( − bω( ade−cd/b b )+cd db , DNc D+lA∗2 ) exists. Proof: First, we will prove that A∗2,N ∗ 2 ∈ R. It is known that all parameters are positive, then a = mlDNce σz, b = lmr2, and d = kz are positive. So, it yields: ade−cd/b b > 0 > − 1 e . Let y = ade −cd/b b . According to the nature of the Lambert ω function, since y > −1 e , then ω(y) ∈ R. So that A∗2,N ∗ 2 ∈ R. Next, we will prove A∗2,N ∗ 2 ≥ 0. It is known that a ≤−c, then because b,d > 0, it is obtained as follows: a ≤−c ⇒ ade−cd/b b ≤− cde−cd/b b . Let x = −cd b , then y = ade x b . Let g = xex, as follows: g = xex = − cde−cd/b b . Therefore, y ≤ g. Also, we have x = ω(g). Since 0 < a ≤ −c, then y and g are ascending (increas- ing) functions. Then we apply Lambert ω function as follows: ω(y) ≤ ω(g) ⇒ ω (ade−cd/b b ) ≤− cd b ⇒− bω ( ade−cd/b b ) + cd b ≥ 0 ⇒−d bω ( ade−cd/b b ) + cd bd ≥ 0 ⇒ dA∗2 ≥ 0. Since d > 0, then A∗2 ≥ 0. Furthermore, because all parameters are positive, then N∗2 = NcD D+lA∗2 ≥ 0. Therefore, the second equilibrium point exists. IV. STABILITY ANALYSIS In this section, we will discuss the stability of the equilibrium points that had been obtained. Further on we use A = A(t) and N = N(t) for brevity. Define dA/dt = g1(A,N) and dN/dt = g2(A,N). Note the following: ∂ ∂A g1(A,N) = −Ar1kzlNe(σ+kA)z (lNe(σ+kA)z + r2)2 + r1N lNe(σ+kA)z + r2 −m ∂ ∂N g1(A,N) = Ar1r2 (lNe(σ+kA)z + r2)2 ∂ ∂A g2(A,N) = −lN ∂ ∂N g2(A,N) = −D − lA. By using the above equations, the Jacobian matrix is obtained as follows, J(A,N) =   ∂∂Ag1(A,N) ∂∂N g1(A,N) ∂ ∂A g2(A,N) ∂ ∂N g2(A,N)   . Theorem 5. Given all of the parameters are positive. If ∂g1(A∗1,N ∗ 1 )/∂A < 0 ⇒ r1Nc lNceσz+r2 < m, then the first equilibrium point (A∗1,N ∗ 1 ) is stable. Proof: Note that ∂ ∂N g1(A ∗ 1,N ∗ 1 ) = 0. Substituting the first equilibrium point into the Jacobian matrix, it yields: J(A∗1,N ∗ 1 ) =   ∂∂Ag1(A∗1,N∗1 ) 0 ∂ ∂A g2(A ∗ 1,N ∗ 1 ) ∂ ∂N g2(A ∗ 1,N ∗ 1 )   . According to the Jacobian matrix, the polynomial char- acteristic is as follows: |J(A∗1,N ∗ 1 ) −λI| = 0( ∂ ∂A g1(A ∗ 1,N ∗ 1 ) −λ )( ∂ ∂N g2(A ∗ 1,N ∗ 1 ) −λ ) = 0. Furthermore, the eigenvalues of the Jacobian matrix at the first equilibrium point are obtained: λ∗1 = ∂ ∂A g1(A ∗ 1,N ∗ 1 ) = r1Nc lNceσz + r2 −m λ∗2 = ∂ ∂N g2(A ∗ 1,N ∗ 1 ) = −D. Because all parameters are positive and r1Nc lNceσz+r2 < m ⇒ r1Nc lNceσz+r2 − m < 0, then λ∗1 < 0 and λ∗2 < 0. So, the first equilibrium point is stable. Theorem 6. Given all of the parameters are positive. If ∂g1(A∗2,N ∗ 2 )/∂A < 0, then the second (A ∗ 2,N ∗ 2 ) equilibrium point is stable. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 5/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight Proof: By substituting the second equilibrium point into the Jacobian matrix, it yields: J(A∗2,N ∗ 2 ) =   ∂∂Ag1(A∗2,N∗2 ) ∂∂N g1(A∗2,N∗2 ) ∂ ∂A g2(A ∗ 2,N ∗ 2 ) ∂ ∂N g2(A ∗ 2,N ∗ 2 )   . Let: J(A∗2,N ∗ 2 ) =  t1 t2 t3 t4   . Furthermore, the polynomial characteristic is |J(A∗2,N ∗ 2 ) −λI| = 0, λ2 − (t1 + t4)λ + t1t4 − t2t3 = 0. So, the eigenvalues are obtained: λ∗∗1,2 = t1 + t4 ± √ (t1 + t4)2 − 4(t1t4 − t2t3) 2 . Note that all parameters are positive. Then: t1 < 0, t2 > 0, t3 < 0, t4 < 0, since by the theorem we have t1 < 0. Therefore the real part is negative: Re(λ∗∗1,2) < 0. So, the second equilibrium point is stable. V. BIFURCATION ANALYSIS In the previous section, the equilibrium point stability conditions have been obtained. In this section, we will discuss bifurcation analysis of the model (1) and (2). Since the equilibrium point stability conditions are a function of parameters, then these functions divide the parameter space into several regions. Furthermore, the stability areas are formed on the parameter space. First, we investigate the bifurcation for the varied Nc and D parameters. To find the bifurcation line, at least one of the eigenvalues in Theorem 5 should be 0. Let: λ∗1 = ∂ ∂A g1(A ∗ 1,N ∗ 1 ) = r1Nc lNceσz + r2 −m = 0, λ∗2 = ∂ ∂N g2(A ∗ 1,N ∗ 1 ) = −D 6= 0. It means that a = −c. Hence, we get the first equi- librium point A∗1 = 0 and N ∗ 1 = Nc. Since the other parameter values have been determined in Table I, then we get the bifurcation line as function of parameters: B(Nc,D) = r1Nc lNceσz + r2 −m. The graph of the function B(Nc,D) = 0 can be seen in Figure 1. This function divides the parameter space into two regions, i.e., D1 and D2. In the D1 region for (Nc,D) = (0.048, 1) and D2 region for (Nc,D) = (0.378266, 1), the phase portrait can be seen in Figure 2 and Figure 3 respectively. In both figures, the red dot is the first equilibrium point (A∗1,N ∗ 1 ), while the blue dot is the second equilibrium point (A∗2,N ∗ 2 ). In the region D1 for the parameter value Nc = 0.048 and D = 1, the solution A(t) converges at 0 while the solution N(t) converges at 0.048. So, the solution (A(t),N(t)) is stable at the first equilibrium point (A∗1,N ∗ 1 ) = (0, 0.048), as can be seen in Figure 2. Note that in the phase portrait, the second equilibrium point (A∗2,N ∗ 2 ) = (−1.865, 0.077) is not only non-existent but also unstable. So, the microalgae population tends to be 0 g cm−3 and nutrition tends to be 0.048 g cm−3 as t →∞. In the region D2 for the parameter value Nc = 0.3783 and D = 1, the solution A(t) converges at 9.1807 while the solution N(t) converges at 0.1334. So, the solution (A(t),N(t)) is stable at the second equilibrium point (A∗2,N ∗ 2 ) = (9.1807, 0.1334), as can be seen in Figure 3. Note that in the phase portrait, the first equilibrium point (A∗1,N ∗ 1 ) = (0, 0.3783) is unstable. So, the microalgae population tends to be 9.1807 g cm−3 and nutrition tends to be 0.1334 g cm−3 as t →∞. Based on the bifurcation analysis for the varied Nc and D parameters, it shows that there is a shift in the stability of the equilibrium point. In the region D1, the first equilibrium point is stable while the second equi- librium point is unstable. In the region D2, the second equilibrium point is stable while the first equilibrium point is unstable, thus a transcritical bifurcation occurs. In the next discussion, we investigate the bifurcation for the varied I0 and z parameters. To find the bifurca- tion line, at least one of the eigenvalues in Theorem 6 should be 0. Let: λ∗∗1 = t1 + t4 + √ (t1 + t4)2 − 4(t1t4 − t2t3) 2 6= 0, λ∗∗2 = t1 + t4 − √ (t1 + t4)2 − 4(t1t4 − t2t3) 2 = 0. We easily derive: λ∗∗2 = t1t4 − t2t3 = AlNr1r2 + A(D + lA)r1kzlNe (σ+kA)z (lNe(σ+kA)z + r2)2 − r1N(D + lA) lNe(σ+kA)z + r2 + m(D + lA) = 0. It means that a < −c. Hence, we get the second equi- librium point (A∗2,N ∗ 2 ) defined earlier in Section III. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 6/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight Table I: Parameter table with variations Nc and D. Parameter Description Value Reference D Dilution (0,1]day−1 assumption m Microalgae natural death rate 0.3day−1 [8] k Algae turbidity 0.032cm−1 [5] β Proportional constant 10J−1 cm3 [5] σ Water turbidity 0.1cm−1 [5] z Chemostat depth 8cm assumption l Nutrition uptake rate 0.2day−1 [5] I0 Sunlight intensity 50J cm−1 day−1 [5] Nc Stock nutrition (0,1]g cm−3 assumption q Minimum quota cell for algae grow 0.05g cm−3 [5] Table II: Parameter table with variations I0 and z. Parameter Definition Value Reference D Dilution 0.24day−1 [5] m Microalgae natural death rate 0.3day−1 [8] k Algae turbidity 0.032cm−1 [5] β Proportional constant 10J−1 cm3 [5] σ Water turbidity 0.1cm−1 [5] z Chemostat depth (0,100]cm assumption l Nutrition uptake rate 0.2day−1 [5] I0 Sunlight intensity (0,100]J cm−1 day−1 assumption Nc Stock nutrition 1g cm−3 assumption q Minimum quota cell for algae grow 0.05g cm−3 [5] Since the other parameter values have been determined in Table II, then we get the bifurcation line as function of parameters: B(I0,z) = λ ∗∗ 2 . The graph of the function B(I0,z) = 0 can be seen in Figure 4. This function divides the parameter space into two regions, i.e., D1 and D2. In the D1 region for (I0,z) = (88, 50) and D2 region for (I0,z) = (88, 20), the portrait phase can be seen in Figure 5 and Figure 6 respectively. In both figures, the red dot is the first equilibrium point (A∗1,N ∗ 1 ) while the blue dot is the second equilibrium point (A∗2,N ∗ 2 ). In the region D1 for the parameter value I0 = 88 and z = 50, the solution A(t) converges at 0 while the solution N(t) converges at 1. So, the solution (A(t),N(t)) is stable at the first equilibrium point (A∗1,N ∗ 1 ) = (0, 1), as can be seen in Figure 5. Note that in the phase portrait, the second equilibrium point (A∗2,N ∗ 2 ) = (−0.3216, 1.3661) is not only non-existent but also unstable. So, the microalgae population tends to be 0 g cm−3 and nutrition tends to be 1 g cm−3 as t →∞. In the region D2 for the parameter value I0 = 88 and z = 20, the solution A(t) converges at 3.437 while the solution N(t) converges at 0.259. So, the solution (A(t),N(t)) is stable at the second equilibrium point (A∗2,N ∗ 2 ) = (3.437, 0.258), as can be seen in Figure 6. Note that in the phase portrait, the first equilibrium point (A∗1,N ∗ 1 ) = (0, 1) is unstable. So, the microalgae population tends to be 3.437 g cm−3 and nutrition tends to be 0.258 g cm−3 as t →∞. Based on the bifurcation analysis for the varied I0 and z parameters, it shows that there is a shift in the stability of the equilibrium point. In the D1 area, the first equilibrium point is stable while the second equi- librium point is unstable. For the D2 region, the second equilibrium point is stable while the first equilibrium point is unstable, thus a transcritical bifurcation occurs. VI. CONCLUSION According to the bifurcation analysis in the previous section, the transcritical bifurcation occurred for the varied Nc and D parameters. In the same way, the transcritical bifurcation occurred for the varied I0 and z parameters. Furthermore, for the varied Nc and D parameters, the occurrence of transcritical bifurcations indicated that microalgae have the minimum supply of nutrition to survive. In the other words, microalgae can survive, if nutrition supply is more than the minimum Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 7/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight 0.0 0.2 0.4 0.6 0.8 Nc 0.0 0.2 0.4 0.6 0.8 D D1 D2 Parameter Diagram Fig. 1: Parameter Diagram (Nc,D). −2.0 −1.5 −1.0 −0.5 0.0 A −0.4 −0.2 0.0 0.2 0.4 N Phase Portrait Fig. 2: Phase portrait in region D1 for Nc = 0.048, D = 1. 0 2 4 6 8 10 A −0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 N Phase Portrait Fig. 3: Phase portrait in region D2 for Nc = 0.3783, D = 1. 0 10 20 30 40 50 60 70 I0 0 10 20 30 40 50 60 70 z D2 D1 Parameter Diagram Fig. 4: Parameter Diagram (I0,z). −0.4 −0.2 0.0 0.2 0.4 A 0.8 1.0 1.2 1.4 1.6 N Phase Portrait Fig. 5: Phase portrait in region D1 for I0 = 88, z = 50. −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 A 0.0 0.2 0.4 0.6 0.8 1.0 N Phase Portrait Fig. 6: Phase portrait in region D2 for I0 = 88, z = 20. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 8/9 https://doi.org/10.55630/j.biomath.2023.01.307 Mahardhika et al., Bifurcation analysis of mathematical model of microalgae growth under the influence of sunlight supply of nutrition required for their survival. In this study, microalgae can be sustained, if the supply of nutrition is greater than 0.0967 g cm−3. Meanwhile, for the varied I0 and z parameter, the occurrence of transcritical bifurcation can be interpreted as follows. Microalgae that live in chemostat at depth z requires sufficient light intensity to survive. In other words, microalgae can be sustained in the chemostat at a constant light intensity I0, if the chemostat has a depth no more than a certain value. Thus, microalgae have the minimum light intensity required for their survival in a chemostat at a certain depth. VII. ACKNOWLEDGEMENTS The authors would like to thank the Mathematics Department of Gadjah Mada University and Mr. Hardi who have encouraged the authors to write this research paper. REFERENCES [1] C.-Y. Chen, K.-L. Yeh, R. Aisyah, D.-J. Lee, J.-S. Chang, Culti- vation, photobioreactor design and harvesting of microalgae for biodiesel production: a critical review, Bioresource Technology, 102(1):71–81, 2011. [2] P. Basu, Chapter 12: Production of synthetic fuels and chemicals from biomass, Biomass Gasification, Pyrolysis Torrefaction. Practical Design and Theory, Third Edition, Elsevier, 2018. [3] I. A. Severo, M. C. Deprá, L. Q. Zepka, E. Jacob-Lopes, Chapter 8: Carbon dioxide capture and use by microalgae in photobioreactors, Bioenergy with Carbon Capture and Storage, Elsevier, 2019. [4] Md. S. Alam, Md. S. Tanveer, Chapter 5: Conversion of biomass into biofuel: a cutting-edge technology, Bioreactors. Sustainable Design and Industrial Applications in Mitigation of GHG Emissions, Elsevier, 2020. [5] E. O. Alzahrani, M. M. El-Dessoky, P. Dogra, Global dynamics of a cell quota-based model of light-dependent algae growth in a chemostat, Communications in Nonlinear Science and Numerical Simulation, 90:105295, 2020. [6] M. R. Droop, Some thoughts on nutrient limitation in algae, Journal of Phycology, 9(3):264–272, 1973. [7] A. Rodger, Concentration determination using Beer-Lambert law, Encyclopedia of Biophysics, Springer, Berlin, 2018. [8] R. Serra-Maia, O. Bernard, A. Gonçalves, S. Bensalem, F. Lopes, Influence of temperature on Chlorella vulgaris growth and mortality rates in a photobioreactor, Algal Research, 18:352–359, 2016. [9] S. D’Angelo, L. Gabrielli, L. Turchet, Fast approximation of the Lambert W function for virtual analog modelling, Proceedings of the 22nd International Conference on Digital Audio Effects, Birmingham, UK, 2019. [10] L. Perko, Differential Equations and Dynamical Systems, Third Edition, Springer, New York, 1991. [11] G. J. Olsder, J. W. van der Woude, Mathematical Systems Theory, Third Edition, Delft University Press, 2005. [12] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Third Edition, Springer, New York, 2004. [13] G. C. Layek, An Introduction to Dynamical Systems and Chaos, Springer, New Delhi, 2015. [14] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, 2013. [15] C. A. Laamanen, J. A. Scott, Chapter 1: Microalgae bio- fuel bioreactors for mitigation of industrial CO2 emissions, Bioreactors. Sustainable Design and Industrial Applications in Mitigation of GHG Emissions, Elsevier, 2020. [16] M. A. Vale, A. Ferreira, J. C. M. Pires, A. L. Gonçalves, Chapter 17: CO2 capture using microalgae, Advances in Carbon Capture. Methods, Technologies and Applications, Elsevier, 2020. [17] H. A. Frank, R. J. Cogdell, Chapter 8.6: Light Capture in Photosynthesis, Comprehensive Biophysics, 8, Elsevier, 2012. [18] B. Henderson-Sellers, H. R. Markland, Decaying Lakes: The Origins and Control of Cultural Eutrophication, John Wiley and Sons, New York, 1987. [19] J. M. Hollas, Modern spectroscopy, John Wiley and Sons, New Jersey, 2004. [20] R. G. de Vega, S. Goyen, T. E. Lockwood, P. A. Doble, E. F. Camp, D. Clases, Characterisation of microplastics and unicellular algae in seawater by targeting carbon via single particle and single cell ICP-MS, Analytica Chimica Acta, 1174:338737, 2021. [21] T. P. Dence, A brief look into the Lambert W function, Applied Mathematics, 4(6):887–892, 2013. Biomath 12 (2023), 2301307, https://doi.org/10.55630/j.biomath.2023.01.307 9/9 https://doi.org/10.1016/j.biortech.2010.06.159 https://doi.org/10.1016/j.biortech.2010.06.159 https://doi.org/10.1016/j.biortech.2010.06.159 https://doi.org/10.1016/j.biortech.2010.06.159 https://www.elsevier.com/books/biomass-gasification-pyrolysis-and-torrefaction/basu/978-0-12-812992-0 https://www.elsevier.com/books/biomass-gasification-pyrolysis-and-torrefaction/basu/978-0-12-812992-0 https://www.elsevier.com/books/biomass-gasification-pyrolysis-and-torrefaction/basu/978-0-12-812992-0 https://doi.org/10.1016/B978-0-12-816229-3.00008-9 https://doi.org/10.1016/B978-0-12-816229-3.00008-9 https://doi.org/10.1016/B978-0-12-816229-3.00008-9 https://doi.org/10.1016/B978-0-12-816229-3.00008-9 https://doi.org/10.1016/B978-0-12-821264-6.00005-X https://doi.org/10.1016/B978-0-12-821264-6.00005-X https://doi.org/10.1016/B978-0-12-821264-6.00005-X https://doi.org/10.1016/B978-0-12-821264-6.00005-X https://doi.org/10.1016/j.cnsns.2020.105295 https://doi.org/10.1016/j.cnsns.2020.105295 https://doi.org/10.1016/j.cnsns.2020.105295 https://doi.org/10.1016/j.cnsns.2020.105295 https://doi.org/10.1111/j.1529-8817.1973.tb04092.x https://doi.org/10.1111/j.1529-8817.1973.tb04092.x https://doi.org/10.1007/978-3-642-35943-9_775-1 https://doi.org/10.1007/978-3-642-35943-9_775-1 https://doi.org/10.1016/j.algal.2016.06.016 https://doi.org/10.1016/j.algal.2016.06.016 https://doi.org/10.1016/j.algal.2016.06.016 https://doi.org/10.1016/j.algal.2016.06.016 http://www.lucaturchet.it/PUBLIC_DOWNLOADS/publications/conferences/Fast_approximation_of_the_lambert_w_function_for_virtual_analog_modelling.pdf http://www.lucaturchet.it/PUBLIC_DOWNLOADS/publications/conferences/Fast_approximation_of_the_lambert_w_function_for_virtual_analog_modelling.pdf http://www.lucaturchet.it/PUBLIC_DOWNLOADS/publications/conferences/Fast_approximation_of_the_lambert_w_function_for_virtual_analog_modelling.pdf http://www.lucaturchet.it/PUBLIC_DOWNLOADS/publications/conferences/Fast_approximation_of_the_lambert_w_function_for_virtual_analog_modelling.pdf https://books.google.com/books?id=VFnSBwAAQBAJ https://books.google.com/books?id=VFnSBwAAQBAJ http://www.boekhandelkrings.nl/images/boeken/90/656/2/2/9789065622136.pdf http://www.boekhandelkrings.nl/images/boeken/90/656/2/2/9789065622136.pdf https://doi.org/10.1007/978-1-4757-3978-7 https://doi.org/10.1007/978-1-4757-3978-7 https://doi.org/10.1007/978-81-322-2556-0 https://doi.org/10.1007/978-81-322-2556-0 https://books.google.com/books?id=XYIpBAAAQBAJ https://books.google.com/books?id=XYIpBAAAQBAJ https://doi.org/10.1016/B978-0-12-821264-6.00001-2 https://doi.org/10.1016/B978-0-12-821264-6.00001-2 https://doi.org/10.1016/B978-0-12-821264-6.00001-2 https://doi.org/10.1016/B978-0-12-821264-6.00001-2 https://doi.org/10.1016/B978-0-12-819657-1.00017-7 https://doi.org/10.1016/B978-0-12-819657-1.00017-7 https://doi.org/10.1016/B978-0-12-819657-1.00017-7 https://doi.org/10.1016/B978-0-12-819657-1.00017-7 https://doi.org/10.1016/B978-0-12-374920-8.00808-0 https://doi.org/10.1016/B978-0-12-374920-8.00808-0 https://doi.org/10.2134/jeq1989.00472425001800010030x https://doi.org/10.2134/jeq1989.00472425001800010030x https://doi.org/10.2134/jeq1989.00472425001800010030x https://books.google.com/books?id=lVyXQZkcKKkC https://books.google.com/books?id=lVyXQZkcKKkC https://doi.org/10.1016/j.aca.2021.338737 https://doi.org/10.1016/j.aca.2021.338737 https://doi.org/10.1016/j.aca.2021.338737 https://doi.org/10.1016/j.aca.2021.338737 https://doi.org/10.1016/j.aca.2021.338737 https://doi.org/10.4236/am.2013.46122 https://doi.org/10.4236/am.2013.46122 https://doi.org/10.55630/j.biomath.2023.01.307 Introduction Microalgae Growth Model Microalgae Microalgae growth models under influence of sunlight Determine Equilibrium Points Stability Analysis Bifurcation Analysis Conclusion Acknowledgements References