Abstract
Though there have been many attempts to address growth kinetics in algal photobioreactors, surprisingly little have attempted an agent-based modelling (ABM) approach. ABM has been heralded as a method of practical scientific inquiry into systems of a complex nature and has been applied liberally in a range of disciplines including ecology, physics, social science, and microbiology with special emphasis on pathogenic bacterial growth. We bring together agent-based simulation with the Photosynthetic Factory (PSF) model, as well as certain key bioreactor characteristics in a visual 3D, parallel computing fashion. Despite being at small scale, the simulation gives excellent visual cues on the dynamics of such a reactor, and we further investigate the model in a variety of ways. Our parallel implementation on graphical processing units of the simulation provides key advantages, which we also briefly discuss. We also provide some performance data, along with particular effort in visualisation, using volumetric and isosurface rendering.
1. Introduction
The motivation for optimising growth kinetics in phytoplankton, specifically algae, is rooted in their use for a variety of purposes including dietary supplements such as spirulina [1], and Astaxanthin [2]. Astaxanthin in particular, is a valuable carotenoid often used for pigmentation in salmon and trout, as well as for human consumption, due to its antioxidant qualities [3]. The mass cultivation of algae is not only done for extracting dietary supplements. It is also performed for other tasks such as effluent treatment [4] and biodiesel production (though this is still in infancy) [5].
Previous agent-based algal growth models typically assume a 1-dimensional lattice [6, 7]. Our work focuses on a 2-dimensional lattice in an attempt to better model the local interactions of hydrodynamics in order to more accurately determine illumination history. Illumination history is important to the cell division rates in a culture [8–10], and factors such as fluid dynamics determine the effects of mutual shading between cells, as well as their exposure (or overexposure) to the illumination source (in our work, we assume a parallel light source from either side of the lattice). In our simulation, this combination of photolimitation, photoinhibition, and mutual shading determines the illumination history of a cell. We attempt to capture these factors with some accuracy in counting the state transitions between specific states in the PSF model.
In our previous work [11, 12], we have considered preliminary modelling of semirealistic photo penetration into a 2-dimensional flat-panel photobioreactor, using a reinterpreted model of ferromagnetism devised by Kawasaki [13] with a growth mechanism. The second modified model we proposed [11] improved upon the first [12] by including the work of Eilers and Peeters towards a Markov process describing the photosynthetic cycle of phytoplankton [14–16]. We use this Photosynthetic Factory (PSF) model to decide cell division rates per cell. We provide more details on the PSF model in Section 2.
At this point, we have established motivation for cultivating algae at large scale, and some previous work towards modelling it. In our work, we restrict ourselves only with rectilinear, bubble-column photobioreactors at very small scale. This is by no means the only cultivation method; however, there are also dome-shaped reactors, flat panel airlift reactors, and tubular reactors [17]. Our reason for choosing the bubble-column reactor (also known as a vertical-column bioreactor) to simulate is pragmatic. These reactors scale reasonably well from laboratory scale to indoor mass production scale [18]. The emphasis on scaling in our choice of simulated bioreactor is also to consider the prohibitive cell counts that are typically encountered in even small-scale bioreactors. A lattice capable of accommodating a cellular density of mL−1 for a good yield [2] would be highly CPU and memory intensive and that would be simply for one millilitre.
Using computational methods to simulate growth is not a new concept. Simulation packages such as BSim [19] and BacSim [20] indicate a maturity in the area, though these two packages are primarily intended for pathogenic growth simulation and visualisation. In the past, quantitative analyses of the dynamics of growth in microorganisms were predominantly based on models using ordinary differential equations [9]. While adequate, these models do not consider the system on a microlevel. They instead focus on macrolevel, where the input and output of a system in terms of growth rate and other factors match with a set of differential equations, to a certain degree.
Figures 1 and 2 show examples of our simulation, demonstrating, the complex patterns that emerge from local interaction.
Cellular automata and agent-based methods have proven themselves to be a good method for investigating less mathematically tractable problems [21], for example, the ubiquitous simulation of flocking behaviour in birds and fish [22]. While differential equations were not intuitive in analysing such behaviour, the agent-based computational model is trivially simple. It is so elegant; it has been used very successfully in films all over the world [23]. In essence, the ABM methodology is generative: the macrolevel phenomenon is induced from the collective behaviour of interacting autonomous agents. There are some difficulty in taking this approach, since describing the system on a microlevel demands a thorough understanding of the system, from which a set of behaviours can be inferred. A typical modelling exercise in this light involves identifying the agents and rules and then determining the environment and communication topologies [21].
Not many research models have been proposed for growth in microalgae systems. We can only speculate that this is due to the prohibitive numbers of cells at play [2], which make it simply impractical to simulate. However, it is possible to simulate these at a cluster level rather than individual cell level. Scheffer et al. proposed in 1995 the concept of a “superindividual” for exactly this reason [24]. While still an approximation, it gives some insight into the possible dynamics of prohibitively large systems. Our work presents a spearheaded effort to simulate these systems by exploiting the computing power of graphical processing units (GPUs) [25] and hopefully to circumvent the use of such approximations. We hypothesise there to be a synergetic advantage in cellular level individualistic models in this domain. We attempt to create a model capable of roughly emulating a small-scale (laboratory bench scale) agent-based simulation of algal growth, in order to shed light on growth kinetics and spatial patterns that may emerge and give rise to the growth dynamics other models often assume.
GPUs are not always intuitive to utilise in a simulation. In recent years, NVIDIA has released a framework named Compute Unified Device Architecture (CUDA) making their use in general-purpose programming much more streamlined [26]. In the past, fragment shaders were used to harness the GPU for general-purpose use [27]. This practice and use of the newer CUDA framework has become known as general-purpose GPU (GPGPU). In our work, we use NVIDIA’s CUDA for accelerating simulations.
Much of the efforts expended toward accommodating GPUs are spent satisfying hardware requirements. CUDA devices have a hierarchy of memory banks and caches, with varying access penalties and scope. In addition, the architecture imposes several restrictions on the configuration of threads, and how device-specific code is organised into “kernels,” resembling C functions with some special syntax.
In Section 2, our method is presented along with a brief overview of the constituents of the algorithm and how this builds upon our earlier work. We present results from a few experiments in Section 3, along with visualisation methods. Finally, we discuss and conclude in Sections 4 and 5.
2. Method
A thorough discussion of the PSF model is out the scope of this paper. However, we provide a brief summary for self-containment. The Photosynthetic Factory model facilitates a simplification of the photosynthesis process by considering it as a Markov process [14]. Three discrete states are allowed: activated, resting, and inhibited. Cells may transition between these states with certain predefined probabilities. The state diagram is shown in Figure 3. Our modifications are simply the addition of cell division and increment states. The former allows the cell to split depending on , and the latter increments a counter. The counter is used to influence the growth rate, and is described later.
Our method draws from our previous work and expands upon several aspects. Firstly, we believe that geometry is important to a certain degree (especially in the light of hydrodynamics playing a part in illumination history of cells), and as such, we have adopted a three-dimensional lattice. Secondly, the PSF model does not account for negative biomass growth, which we have added somewhat arbitrarily as a low cell death probability in case of low and excessive illumination. We also ensure that cell division rates are correlated approximately linear to the number of state transitions from the active to resting states. This is due to supporting literature suggesting that there is a correlation here [7]. We have done this using a logarithmic function to ensure a valid probability:
Values of above are typically in the range of 0 to 200 in our simulations.
Another improvement upon our previous model is our photopenetration characteristics. Photopenetration into a bioreactor vessel is not a trivial consideration. Irradiance on cells is not constant throughout operation of the reactor. Hydrodynamics of bubbles and gravity through the medium cause cells to move and sporadically shadow or being shadowed by other cells to a certain extent. We modify light penetration by computing a light map due to mutual shading in the reactor using CUDA. This is then used by the cells to determine the likely irradiance available to them in a separate set of CUDA kernels.
Apart from the changes directly to the model itself, we also add bubbling effects, which are prevalent in bubble column reactors. Our method for accomplishing this is to randomly release a third species of site “spin” near the bottom of the vessel, which randomly diffuse upwards. This is to include the effects of bubbles rising from the sparger in a real photobioreactor. These bubbles act as an upward, monopoled repulsive particle to other species in the tank. Mass transfer is mostly studied for reducing how long it takes a vessel to reach a homogenous mixture of nutrients [18]. At this point, we are not considering the transfer of nutrients, however.
Using a reinterpreted ferromagnetism model for a growth simulation may be somewhat unorthodox, but in our previous work, this has resulted in a fairly realistic clustering of algal cells and random diffusion in the vessel. Coupled with a weak gravity force [28], and the lattice-gas characteristics of the Kawasaki exchange model, we have no reason to disregard this approach to a solution. Given this gravity force, it has been previously shown that the Kawasaki exchange model is effective in modelling sedimentation [28] specifically for species of different mass.
We endeavour to keep assumptions to a minimum, but unfortunately some must be made. We assume nutrient saturation in the vessel at all times in order to ensure that irradiance is the only variable, with illumination history resulting from stochastic hydrodynamic diffusion. The reason that illumination is frequently studied in isolation is due to its overarching influence in growth kinetics, and especially so when the culture is not under any other significant stress [8, 29] (NaCL, heat, and nutrient deficiency). It is for this reason that we consider only the influence of photopenetration through the medium in the cultivation vessel. Other simplifying assumptions are also made, including the dynamics of bubbles. This is somewhat simplified since we do not allow nucleation of bubbles in the lattice.
The probability of cell division is key to the simulation and depends directly on the number of state transitions a cell has had between activated and resting states. Illumination is also transitively important: the number of state transitions between activated and resting depends on how much illumination the cell receives.
For modelling light penetration into the vessel, we use a standard exponential absorption function in the spirit of the Lambert-Beer law from the study of optics [30, 31]. Scatterance is another aspect to hydrological optics, in tandem with photoabsorption [32]; however, in this work, we focus on absorption in the medium by a constant attenuance exponent () and manually compute mutual cell shading.
The full equation describing illumination across the bioreactor vessel is shown in
Here, cell shading is represented by and . These are simply the numbers of cells between the current cell and the left and right vessel walls. is the position of the cell and is the width of the reactor. and provide attenuation and scale to the degradation of the light penetration.
Finally, the actual formula for is simply a scaled from (1).
Figure 4 gives a visual representation of illumination from the left to the right sides of the vessel. represents the left edge, and the right edge. While this surface describes the general availability of irradiance on cells, it does not take into account the effects of mutual shading. Mutual shading applies a further fractional coefficient to the illumination, and its landscape is shown in Figure 5. The “shading” axis represents how much of the vessel is obstructed for the cell under consideration. The inflection at divides the left edge from the right edge.
General pseudocode for the parallel simulation is shown in Algorithm 1. Race conditions are especially difficult to deal with in a parallel-lattice-based simulation of spin exchanges. We have therefore opted to divide the entire lattice in cubes, each of which is evaluated by one thread, times synchronously. We omit an extensive discussion on the parallel implementation. More details on this are available elsewhere [11]. The effects of mutual shading are computed by executing a CUDA kernel with a , grid of threads, sweeping over the lattice from to .
| ||||||||||||||||||||||||||||||||||||||||||
Algorithms 2 and 3 represent the implementation of the Markov process and the details of the Kawasaki site exchange computations, respectively. The initialisation shown at the start of Algorithm 1 involves “inoculating” the reactor with 10 randomly positioned cells in resting state with zero-transition counts.
| ||||||||||||||||||||||||
| ||||||||||||||||||||||
The Metropolis probability depends on a temperature variable, which we have fixed to unless otherwise noted. Other parameters used in the simulation are shown in Table 1. Our choice of the PSF parameters , , , and is somewhat arbitrary and not fully appropriate in terms of the original PSF model. The novelty in our approach lies in the use of a pseudo-PSF model in order to extract a spatially heterogeneous, real-time growth probability. We demonstrate this as an agent-based model and accelerate this using parallel computing.
It is prudent to point out that it is unrealistic to assume that all agent-based systems can be modelled as elegantly as “boids” [22]. We do concede that models with excessive parameters are limited in their prediction capabilities to a certain extent. Eilers and Peeters also conceded that their updated model of 1993 [16] with six parameters was a large number, presenting possible difficulty in fitting to collected data. We attempt to mitigate this somewhat by holding all parameters constant. There may be scope in the future for using metaheuristics [33] in fitting data. Metaheuristic algorithms are particularly well suited to dealing with problems containing many parameters. At this point, however, we resign ourselves to use hand-tuned parameters in our work.
In terms of experimentation, we have carried out some tests to gather characteristics of this bioreactor simulation. We measured global growth rates by collecting fill fractions by time step with illumination varying at certain time steps. The addition of gas bubbles in the medium is investigated. Following these experiments, we propose modifications to the algorithm and discuss them in some detail.
Our work includes the use of parallelisation across GPUs so as to increase system size and pursue more accurate aggregate behaviour from cultures. We measured performance, and we suggest methods of increasing it further.
3. Results
The first observation we made was while overilluminating the vessel, resulting in photoinhibition of the cells. During overillumination, cells that spontaneously form clusters are able to thrive, due to mutual shading reducing illumination to more acceptable levels. An example of this is shown in Figure 6. The CUDA kernel computing the shading map proved to be prohibitively expensive to compute, at around 3.3 ms.
We have previously mentioned a fixed value for the temperature parameter in the Metropolis probability used for random diffusion. In our empirical observations, the temperature parameter actually has an influence on the survival of the culture. Though while it may not be direct in the sense of culture temperature at this point, this parameter does appear to decide some level of hydrodynamic turbulence. Low values of this temperature seemed to cause a rapid decline. We postulate that this is caused by a decline in the number of sites available for cell division, caused by clustering of cells. Currently, should a chosen site be occupied by another cell, cell division will not take place in this model. We then modified the division code to guarantee division to take place, provided that there is at least one empty adjacent site and the division probability is satisfied. We have compared the outcomes of the system with these two configurations in Figure 8. The vessels in these images were modified to be flat panels, in order to best visualise the effects in a thin layer.
An immediate conclusion for improving the simulation is difficult to draw from this diagram. Monte Carlo lattice gas simulations do not typically force a site update to take place, as is the nature in these simulations. It is clear, however, that there is a large increase in growth rates, created clearly by forcing division to take place. The effects of site shading are not extremely pronounced in these images. The full 3D effect is more difficult to interpret, but for the guaranteed cell division case is shown in Figure 7.
Shown in Figure 9 is the vessel fill ratio with illumination conditions varying at time steps (), (), and (). The purpose of this was to see the effects of underillumination. What is interesting to note is that while there is an immediate decline at , there is a slow increase towards , during which illumination was still prohibitively low. Observations would suggest that this is caused by cell deaths decreasing mutual shading, and hence allowing more light towards the centre of the vessel.
Another observation made is that various values of the light attenuation exponent () resulted in a different spatial configuration of the culture. Lower values of (which correspond to a clear liquid) caused the culture to thrive in the centre of the vessel, due to overillumination at the edges. Higher values of (cloudy liquid) saw the culture adhering to the left and right walls.
We gathered some observations on the effects of bubbles on spatial configurations. These are shown in Figure 10. The image contains visual cues as to how bubbles affect the growth of the culture. A weak gravity force is still present in these simulations, but in the left two images, gravity caused the sparger to become covered in biomass, and therefore unable to produce bubbles. In the second two images, bubbles overcame gravity and the culture formed toward the centre of the tank. It is not clear whether this is a good analog to models of mass transfer in physical bioreactors, for reasons which are discussed in the next section.
Following these results, we modified the simulation so that cell division is also dependent on illumination. This change came following the observation that cells would divide regardless of whether they are currently over- or underilluminated. In a fully grown culture for example, the culture thins in the centre due to photolimitation, and these would previously be immediately replaced by dividing neighbouring cells. This change would also assist in improving the growth characteristics under varying illumination.
It is worth noting at this stage, photolimitation causes cell death in the centre of the reactors, in which case cell division takes place in some circumstances. This is the cause of red (dark) coloured cells near the centre and elsewhere throughout the reactor.
Conveying a sense of qualitative realism is conducive to obtaining more trustworthy quantitative results. Therefore, we have been making a concerted effort in obtaining a recognisable growth pattern from as few a set of parameters and explicit instructions for individual cells. At no point do we specifically control cell division from a centralised location. Cells divide based on their own state traversal in the diagram shown in Figure 3.
Due to the considerable difficulty in interpreting images of the bioreactor, we have experimented with volumetric rendering to improve legibility, using a slightly modified version of the sample in the CUDA SDK which accompanies the CUDA toolkit and programming guide [26]. In Figure 11, we show a comparison of this method with the original. Printed images such as these will invariably be difficult to interpret. However, the simulation itself being interactive is a very useful tool for gathering a good sense of spatial configuration.
(a)
(b)
It is perhaps easier to see the advantage of volumetric rendering in this case by examining Figure 12.
(a)
(b)
The rendering shown in Figure 13 is rendered using isosurface extraction with Phong shading and a marching cubes algorithm implementation in a modified CUDA SDK sample “Marching Cubes” [26]. An appreciation of depth is much easier in this image. Clusters on the edge of the vessel are not closed.
Cursory performance indications were concerning, since for every random site chosen in the sublattice, one CUDA kernel must be executed synchronously across the entire lattice. A total of CUDA kernels (excluding the mutual cell shading map kernel) are computed for one time step. Moreover, each of these individual kernels, while fast for a relatively sparse vessel, becomes more cumbersome as the vessel becomes more dense. While the vessel has a <1% fill ratio, one such kernel executes for a duration of approximately 93 μs, including synchronisation, which accounts for around 2/3 of this. At approximately 80% fill, this increases to around 300 μs, for each of the 27 kernels. This equates to little more than twice the shading map kernel. Given that mutual cell shading is not giving much more realism than an appropriate global illumination attenuance, this result would suggest that computing cell shading manually is simply too expensive to compute manually.
4. Discussion
At this point, we can speculate why ABM has not been applied liberally in the field of algaculture. All of our simulations would be appropriate for studying, at most, 1 mL of undried biomass at 1 : 1 scale. It is possible, however, that such a reactor could be simulated in its entirety across a grid computer, and we will consider a similar model in the future for parallelisation across several GPUs.
Perhaps the greatest downfall of the model we have brought forward is the lack of hydrodynamic diffusion in clustering cells. Clusters appear fixed in space, whereas singular cells diffuse randomly. It is reasonable to expect a cluster to shift in location, except if it has adhered to the wall of the vessel. For a shallow bioreactor simulation, this is less of a concern. Adherence to the reactor walls is more common in bioreactors where relatively little hydrodynamic movement takes place. The temperature parameter also affects exactly what kind of clustering occurs, if any. For this reason, we have sometimes interpreted the temperature parameter of the ferromagnetism model as an indicator of hydrodynamic turbulence.
We have also experimented by adding bubbles in an effort to emulate the effects of hydrodynamics from a gas sparger for nutrient supplementation and mass transfer in the reactor. It may be possible for a sparger to experience difficulty in producing bubbles with biomass adhering to it, although it is unlikely that biomass deposits would reach a stage where bubbles could no longer enter the vessel. It would appear that gravity and bubbles could counteract each other as well.
5. Conclusion
The importance of modelling microbiological phenomena has already been clearly demonstrated, and we have perpetuated the advantages of agent-based modelling to the simulation of small photobioreactors. While this work is somewhat preliminary, we have also experimented with visualisation techniques for gaining better insights into the spatial configurations generated. A combination of the Photosynthetic Factory model, in tandem with a reinterpreted lattice gas algorithm and cell shading, was accelerated using NVIDIA’s CUDA for graphical processing units. The result was visualised in a number of ways, including point rendering, volumetric rendering, and isosurface extraction.
A great deal of the difficulty involved in applying agent-based modelling to algaculture is the excessive numbers of algal cells that need to be modelled. We are in the process of considering a large-scale implementation of the algorithm presented.