Overview
Everything on a mission to space is meticulously calculated. From an exact path trajectory, to fuel weight with a 1% error margin, “precision” is synonymous with space travel. With that precision must come a myriad of complex mathematical equations for a space mission to even have a chance of success against the many unknown variables of space. While this complex engineering is vital to the safety of astronauts, everything else also has to be optimized to peak efficiency for their long term survival. By modeling various systems and determining pitfalls, scientists can predict and optimize systems to provide food and recycle nutrients in a closed loop.
The modeling for this project are made up of four sections, each relating to different parts of the project:
- Clinostat
- Acetate Overflow Path
- Bioreactor
- PHBV Properties
The Clinostat calculations inform Hardware what the ideal RPM ratios are for simulating microgravity and Martian gravity. Acetate Overflow Path describes the rate of acetate production and concentration in module 2, and informs Wet Lab for some ideal experimental parameters. Bioreactor is a simulation of the entire project with an autotrophic acetogen feeding the acetate precursor into a PHB producer. This is to estimate how much PHB can be produced and how long the bioreactor would be able to run for, informing Policy and Practice about production times and labor cost. Lastly, PHBV Properties creates linear regression equations between mol% 3HV and tensile strength, elongation, and elastic modulus. These equations are so the properties of any given bioplastic can be predicted given the percent composition of 3HV and vice versa.
Clinostat
Introduction
Amongst the four fundamental forces in the laws of physics, gravitational force is the least understood. Not only is it incompatible with our current understanding of quantum theory, it is the weakest force. Due to the magnitude of its weakness, mathematical descriptions of gravity were derived from astronomical scales. Other than simulated micro- and hypergravity (such as free falling and centrifugal force), Earth does not offer many natural forms of varied gravity.
The best microgravity research is done in true microgravity conditions, but mathematical models are excellent tools to gain a predictive sense of what will happen before the research is actually conducted. In our project, we chose to use a dual axis clinostat, Spinova, as our microgravity simulator because it offers a more practical way to study the effects of non-Earth gravity on various parts of the PHAntom project. The rotation cancels out some magnitude of the gravitational vector, allowing us to simulate non-earth gravity conditions. By calculating the extent to which the gravitational vector is canceled out, the ideal rotation ratio for a given gravitational acceleration can be determined.
Methods
The mathematical model can be derived from applying a rotation matrix to describe the effect of a frame rotating another frame (Kim, 2016). The full derivation can be found at the, under the “2 axis Derivation” section.
Equation 3
Rotation matrices are established matrices that describe the motion of a reference frame after a rotation (Equation 3). The rotational matrix provides the means of describing the motion of both rotating frames. If only a single frame were rotating, it would be simple to describe its motion using the angular velocity. In this case, the rotation of Spinova’s outer frame is translated to the inner frame. Therefore, we need to apply the rotational matrix to account for the inner frame’s plane of reference being rotated along a perpendicular axis.
Following matrix multiplication, the 2-axis clinostat is found to be described by the following differential equations (Kim, 2016):
Equation 10
The full derivation can be found below or in the referenced paper.
Each differential equation describes the rate of change of an arbitrary vector in the clinostat system. Because we exist in a 3-dimensional system, we need three vectors in order to completely describe the system. These vectors maintain fixed angles relative to the gravitational vector. When gravity's direction changes, they rotate accordingly while preserving their angular relationships. The specific angles don't matter because they remain constant and cancel out in the equations.
The differential equations are coupled, making them difficult and time-consuming to integrate and find an exact analytical solution. It is much easier to compute a numerical solution using Python. The idea is to approximate very small sections of the solution as linear. The rate of change (differential equations) is known, and we can assume that at least for an infinitesimally tiny amount of time, the rate of change is constant. Taking these assumptions, we can predict where the vector will be pointing after that tiny amount of time. By repeating this process for sequential tiny amounts of time, the computer can generate the full path of the vector for a given amount of time. The Python code for this is available on the Software Repository .
After generating the paths of the normal vectors, we need a quantitative way to evaluate which ones resemble microgravity as closely as possible. This can be done by evaluating the time-averaged acceleration, which is mathematically represented by the magnitude of the average vector. The magnitude of the normal vector was initially normalized to 1, the magnitude of the normal vector at any given time should be equal to 1. Gravitational acceleration on earth is about 9.8 m/s, so the vector at any time point would be equal to 1g. In order to take the direction of the vector into consideration, the three individual components of the vector need to be averaged first. This gives the coordinates of the average vector. The magnitude of this average vector will then be the acceleration that the sample experiences, in g’s.
Equation 20
Assumptions
Given that the solution for the governing differential equations is numerical and computer generated, there are a number of assumptions necessary in order to trust the solution:
- The time step is small enough that the error does not significantly compound
- Computer rounding is insignificant
- The period of the vector path is around 600 seconds
The last bullet point is extremely important due to the method of calculating how many g’s are experienced in the system. Since the clinostat is rotating a sample dependent on constant RPMs, the path of the normal vector must be periodic. Since the differential equations have sinusoidal equations in them, there must be periodicity in the solution (the vector path). That means that after some given time, the path must start to repeat.
This is a problem in averaging the vector. The average must be over an integer number of periods to accurately describe the acceleration experienced. If the simulation runs for a period and a half, then the average will be skewed because of that extra half. This is insignificant in experimental trials because as the time increases, the average will get closer and closer to the real acceleration simply due to how many periods have passed. Computationally, that is not feasible because of how long it would take to run the simulation. Therefore, assuming a certain amount of time will average out close enough to a period is very important for getting an accurate measure of acceleration. Since the period for 0.1 rad/sec would be 600 seconds, this value is used as the baseline period for the simulation.
To test the validity of 600 seconds as the period for all angular velocities, all the differential steps are averaged. If the average of the differential steps is close to 0, 600 seconds is a valid estimate for the period. This is because after an entire period, the symmetry of the path indicates that the differential steps would have canceled out. Heat maps can be generated to show how close each of the RPM ratio averages are to 0.
Figure 1. Vector path integration over 1 second (generated using Python).
Figure 2. Vector path integration over 5 seconds (generated using Python).
Figure 3. Vector path integration over 100 seconds (generated using Python).
Figure 4. Vector path integration over 300 seconds (generated using Python).
Figure 5. Vector path integration over 600 seconds (generated using Python).
By visual confirmation of the heat maps in Figures 1-5, 600 seconds is the optimal choice of runtime. For the vector path integration of 600 seconds, the majority of RPM ratios result in integration values 0.002 and less (Figure 5). As the RPMs increase, the integration value starts to increase. Unfortunately, that would mean that the acceleration calculated for higher value RPM ratios may be less accurate. However, note that the maximum value on the colorbar drops from around 0.40 (Figure 2) to around 0.095 (Figure 5). While higher value RPM ratios may be less reliable, a period of 600 seconds is far better than shorter periods.
Results + Discussion
After the differential equations were derived, it is helpful to see a visualization of the different paths that the normal vector can take.
Figure 6. Normal vector paths of all combinations of RPMs 0.1 to 3.9 in x and y direction.
There are some things to note just by visual inspection. Simple RPM ratios like 1:2, 1:3, etc, have simple paths. Multiples of these ratios yield similar paths. For example, 1.3:3.9 has an almost identical path to 0.1:0.3 because their ratios both simplify to 1:3. However, the higher multiple paths seem to be slightly offset. This is likely due to computer rounding and could be fixed by using a smaller timestep. However, this would require more computing power than is available to the team. Rounding errors in large timesteps compound much faster and can cause deviations from the analytical solution. Knowing this information, the ideal microgravity RPM ratio would have to be the most simplified out of the RPM ratios that share the same simplified ratio.
Figure 7. Graph of 5 second vector average (generated using Python).
Figure 8. Graph of 5 second vector average, normalized to Martian gravity (generated using Python).
Table 1. Minimum and maximum integrations for each of the three types of simulation, with the associated RPM ratio.
The heat map color codes the magnitude of the average vector in each of the RPM ratios (Figure 7, 8). The colorbar is set logarithmically to illustrate the magnitude of differences.
When both RPM values are low in Figure 7, the average vector (acceleration) is also low. As the RPM values increase, the average vector increases. When both of the RPM values in the ratio are the same, the calculated acceleration is very high. As seen in Table 1, the minimum average vector magnitude is seen at (0.2, 0.1). The maximum average vector magnitude is included as a testament to the effect of computer rounding. The magnitude of 4.962 would mean that the clinostat would cause the sample to experience almost 5g of gravity. This may in part be due to the assumption of 600 seconds as a period, since averaging over an incomplete period would skew the value (Figure 5). However, it is much more likely that the time step is too large. As the clinostat moves faster, the simulation can no longer approximate the rate of change accurately. This could skew the values of each component so they are not normalized anymore.
Ideally, the bacteria in our project will be able to fix carbon dioxide and produce PHBV on Mars. This means that more important than simulating microgravity is to simulate Martian gravity. Martian gravity is roughly 0.38g, meaning that to find the RPM ratio to best simulate Martian gravity, the RPM ratio must produce an average vector with a magnitude of about 0.38. In Figure 8, a similar heat map is produced. In this map, the colorbar represents the difference from the average vector magnitude and 0.38. As listed in Table 1, the minimum difference is attributed to the rotation ratio of (3.8, 3.4).
Impact on Hardware
The differential equations of vector path from Kim, 2016 along with the derived method of acceleration quantification has two uses:
- Allowing experiments to be run under different gravitational conditions
- Returning the equivalent number of g’s for a given RPM ratio
If experiments need to be run in microgravity, the quantification system can calculate the ideal RPM ratio. Since the project is designed to be implemented on Mars, the system can also produce the ideal RPM ratio for Martian gravity. Slight adjustments can be made for virtually any quantity of gravity that is less than terrestrial. The simulation can also run the reverse process, where a RPM ratio is fed to the system and the vector path is visualized on a graph. With some slight modifications to the code, it can also be changed to calculate the acceleration in g’s.
Summary
After deriving the differential equations that describe the normal vector movement, it is numerically calculated in Python. The average normal vector is then calculated over a period, which was assumed and validated to be 600 seconds. For microgravity, the average normal vector should be close to 0. The RPM ratio closest to 0 is (0.2, 0.1). The simulation is then adjusted and calculated for Martian gravity (0.38 g’s). The ideal RPM ratio for Martian gravity is found to be (3.8, 3.4).
2-Axis Derivation
The equations explained below are taken from the article referenced (Kim, 2016). A condensed pdf version is also available.
Since the clinostat is rotating, φ (phi) and ω (omega) are going to be used, representing angular position and angular velocity respectively. It is easier to write everything in terms of angular velocity, since RPM is what can be controlled on the clinostat. The following conversion will make it easy to write any φ’s in terms of ω.
Equation 1
Another thing that needs to be defined is the total angular velocity.
Equation 2
The normal ω represents the angular velocity of the outer axis (arbitrarily declared to be the z axis). The ω’ (omega prime) is the inner axis angular velocity (x axis, also arbitrary). It is denoted as prime because the inner frame is going to be changed due to the outer axis by a rotation matrix. They are summed to get ωt, representing the total angular velocity. That will be used later to give the coordinate system the direction.
The rotation matrix is given by the following:
Equation 3
This particular rotation matrix describes the motion around the z axis. Since the inner frame is being rotated around the z axis, the rotation matrix needs to be applied to the inner frame’s movements to make sure it aligns with the outer frame. In fact, 3 vectors can be defined specifically to hold direction (e hat, the hat dictating that the magnitude is 1 so only the direction is important). Then, these vectors can be multiplied by the rotation matrix so the direction due to the rotation can be preserved (e hat prime).
Equation 4
Matrix multiplying the rotation matrix against the direction vectors (and rewriting in terms of ω) yields:
Equation 5
Now that the new directions for the inner axis are established, the directions for the angular velocities can be defined.
Equation 6
Since the magnitude of the angular velocity ω is known, multiplying the direction with its respective direction gives the angular velocity vector. The outer frame ω is rotating around the z axis, so it is being multiplied by the z direction. The inner frame ω’ is rotating around the x axis after being rotated by the z axis, so it needs to be multiplied by the adjusted x direction. Therefore, the new angular velocity vectors can be written as the following matrices:
Equation 7
Plugging these back into Eq. 2, the total angular velocity can be found.
Since the magnitude of the angular velocity ω is known, multiplying the direction with its respective direction gives the angular velocity vector. The outer frame ω is rotating around the z axis, so it is being multiplied by the z direction. The inner frame ω’ is rotating around the x axis after being rotated by the z axis, so it needs to be multiplied by the adjusted x direction. Therefore, the new angular velocity vectors can be written as the following matrices:
Equation 8
Since angular velocity represents change in φ, the angular velocity can be multiplied by the normal vector to yield the change in normal vector. Again the normal vector is some fixed angle away from the gravitational vector. Changing the orientation of the system is equivalent to changing φ, which changes the gravitational vector. Since the normal vector has to be some fixed angle away from the gravitational vector, the normal vector will also change due to the change in φ.
Equation 9
After cross multiplying the two vectors, the following differential equations are derived.
Equation 10
Extra 3 axis info
The derivation for a 3 axis clinostat is done out as a proof of concept, should anyone choose to use it.
The difference between the two axis and three axis clinostats is very simple: the extra axis. That would mean that there is another frame of reference.
Equation 11
Adding another frame as the innermost frame means that the original frames will not have to be changed. This also means that the innermost frame will be rotating around the y axis. Since another frame is being added, there is another aspect of the rotational velocity to be added to the total rotational velocity.
Equation 12
Other than those two additions, the rest follows the same style of derivation as the 2 axis clinostat. The first step would be to find the change in reference frame due to the outermost frame, which has already been calculated in Eq. 4, resulting in Eq. 5. In order to get the directions for the new frame of reference, those equations need to be multiplied by another rotational matrix. Since the (now) middle frame is rotating around the x axis, the following matrix is used.
Equation 13
Multiplying the rotational matrix by the middle frame’s directions yields the innermost frame’s directions.
Equation 14
This yields the following equations, very similar to Eq. 5.
Equation 15
The angular velocities can then be calculated using the newfound directions.
Equation 16
Adding together all of the angular velocities yields the total angular velocity.
Equation 17
This can once again be cross multiplied against the arbitrary acceleration vector.
Equation 18
The cross multiplication then yields the following differential equations, which may be used to describe the motion of the acceleration vector in a three axis clinostat.
Equation 19
- An Experimental and Theoretical Approach to Optimize a Three-Dimensional Clinostat for Life Science Experiments | Microgravity Science and Technology. (n.d.). Retrieved May 21, 2025, from https://link.springer.com/article/10.1007/s12217-016-9529-2
- Clary, J. L., France, C. S., Lind, K., Shi, R., Alexander, J. S., Richards, J. T., Scott, R. S., Wang, J., Lu, X.-H., & Harrison, L. (2022). Development of an inexpensive 3D clinostat and comparison with other microgravity simulators using Mycobacterium marinum. Frontiers in Space Technologies, 3. https://doi.org/10.3389/frspt.2022.1032610
Acetate Overflow Path
Introduction
Acetyl-CoA is a universal metabolic intermediate and the primary precursor for PHB/PHBV bioplastic synthesis. Native producers (e.g., Cupriavidus necator, Cupriavidus basilensis, Haloferax mediterranei) channel acetyl-CoA through the pha pathway—β-ketothiolase (PhaA or BktB) and acetoacetyl-CoA reductase (PhaB)—to generate (R)-3-hydroxybutyryl-CoA (3HB-CoA), which PHA synthase (PhaC) polymerizes into PHB. PHBV incorporates 3-hydroxyvaleryl-CoA (3HV-CoA) units, formed when BktB condenses propionyl-CoA + acetyl-CoA to 3-ketovaleryl-CoA, which is then reduced by PhaB to 3HV-CoA, alongside 3HB-CoA.
However, acetyl-CoA is a transient, high-energy thioester that is neither stable nor exportable. To transfer carbon between modules, we overproduce and secrete acetate in the autotrophic strain, then allow the PHBV producer strain to take it up and regenerate acetyl-CoA. This strategy uses a stable, storable intermediate while preserving tight control over polymer-precursor supply.
Production of acetate not only provides the PHBV-producing strain with a stable carbon feedstock, but acetate is also a useful chemical for fertilizer enrichment (enhancing nutrient uptake) and as a precursor in cellulose acetate fabric manufacturing.
By quantifying acetate production (titer, rate, and yield), we can guide procedures in wet lab, including planning media, aeration, and transfer schedules with appropriate margins. Rather than using a dedicated acetogen, our autotrophic E. coli module is engineered to channel fixed carbon into acetate secretion. Conceptually, we split this “acetogen mimic” into two parts: (1) autotrophy, and (2) acetate overproduction. The second module takes advantage of E. coli’s native acetate overflow pathway, enhanced with a recombinant pta–ackA, to drive export under conditions of glucose excess. Establishing the necessary concentration of glucose the bacteria need for the acetate overflow route would be optimal for carrying out acetate overproduction.
Methods
Assuming that acetate overproduction is related to growth in biomass, a system of coupled differential equations can be produced to mathematically model the environment over time.
Figure 1. Schematic representation of simulation. Each arrow represents a differential equation for how individual components of the simulation affect each other.
Based on the general mass balance equation, mass equations for the three components of the simulation can be produced. The following equations are motivated by Anane et. al.
Equation 1 (General mass balance)
The general mass balance equation is dependent on the initial and current values of the substrate xi. The second term, rX, also accounts for how biomass affects the production of a given substrate. The rate r will depend on the specific substrate and how biomass affects its growth.
Biomass is denoted as X, and the rate of change of biomass is given in the equation below, based on the general mass balance.
Equation 2
The initial biomass concentration is low enough to be taken as negligible. That means the first term is overall negative, which makes sense as a container runs out of space for bacteria to keep multiplying when the concentration of bacteria is very high. The rate μ is the rate at which biomass increases. The second term dictates how biomass increases due to that rate and the current concentration of biomass.
The same is true of the glucose concentration (Eq. 3). Due to the acetate overflow route, there needs to be some initial concentration of glucose so that the bacteria will make acetate. Since glucose is not being replenished from any source, the glucose concentration will be depleted over time based on the rate that E. coli metabolizes glucose. This means that the second term will be negative, showing overall decrease in glucose.
Equation 3
Lastly, the production of acetate is taken to be some rate qsA, which is multiplied by the concentration of biomass (Eq. 4). Similarly to the biomass equation, the initial concentration of acetate is 0. As the acetate concentration increases, it inhibits the production of acetate. That is reflected in the first term being negative. When the concentration of acetate gets too large, the bacteria starts breaking it back down as a source of energy.
Equation 4
The rates that are in the differential equations can be described using the following equations.
From Equation 2, the rate μ can be modeled using non-inhibited Monod-type growth. In the following equation, μ is modeled as a function of uptake rates (q) and yield coefficients (Y). All of the specific coefficient definitions and values can be found in Table 1 at the end of this section. Note that uptake rates also have unique definitions and change based on glucose and acetate concentration. Yield coefficients are fixed values.
Equation 5
Equations 6 to 9 are specific uptake rates for glucose (s). All of the following uptake rates are also derived from Monod-type kinetics. In Equation 6, the general uptake rate of glucose is shown. It will end up being some fraction of the maximum uptake rate, depending on the acetate and glucose concentrations. K are affinity constants to describe how well the organism responds to a certain substrate.
Equation 6
Some of the glucose that is taken up is oxidized for regular cell metabolism. This is represented in Equation 7. The portion that isn’t oxidized goes towards the acetate overflow pathway. Metabolism is dependent on oxygen concentration, which is represented in the fraction as a limiting factor.
Equation 7
The glucose used in the overflow route is expressed in the following equation. It is dependent on the maximum production rate of acetate.
Equation 8
Lastly, the intracellular acetate production rate is described using the following equation, similar to Equation 6 describing glucose uptake rate.
Equation 9
Table 1. Comprehensive list of variables organized by equation, along with the definitions and values (Anane, 2017). The bolded number in the Equations column is the equation that the corresponding variable is defined in.
*The values for F, V, and Si are educated guesses. This is elaborated on in the Assumptions section.
Since the differential equations represent the rate of change of the system, they can be used to calculate the rate of change at a specific time in the system. For small amounts of time, the solution can be approximated as linear and the rate of change a constant. A computer simulation utilizes that to numerically calculate a solution. The simulation is a program written in Python that calculates the concentrations of biomass, acetate, and glucose incrementally. It can be found on the Repository. The initial values of differential equations and rates are calculated using the initial values of biomass, acetate, and glucose concentration.
Assumptions
These differential equations are an estimate of the rates of change of the concentrations of biomass, glucose, and acetate. That being said, there are a couple of assumptions that need to be made for the simulation to be ran:
- Perfect Monod-type growth kinetics
- Negligible change in dissolved oxygen
- Some constants need to be assumed
Monod-type kinetics are a mathematical description of observed growth rates under ideal conditions. The equation does not take into account how the bacteria may react to stress. Especially in a new atmosphere and different gravitational conditions, the bacteria may react negatively.
Since organisms need oxygen to metabolize and live, the dissolved oxygen content in growth media is vital. However, since our project assumes that the bacteria will be grown in some sort of container on Mars, there would probably need to be some sort of bubbler sending a constant stream of CO2 and O2 through the system. This means that the change in dissolved oxygen content would be negligible as it would be constantly replenished.
There are a number of constants that are dependent on wet lab information. The volume of the container and the feed are values that would be adjusted based on the wet lab. Initial glucose concentration can be varied to see the resulting acetate concentration graph. In the meantime, these values are assigned appropriate values for the sake of simulation testing.
Results
Since the feed and volume are not known from the wet lab, they are varied to picture several different wet lab set ups.
The simulation was initially run over 5 hours and 50 hours using the set values from Table 1 (Figure 2, 3).
Figure 2. Acetate concentration (g/L) over time for 5 hours.
Figure 3. Acetate concentration (g/L) over time for 50 hours.
The acetate concentration growth over 5 hours looks exponential while the growth over 50 hours is showing the beginnings of a plateau. This indicates that 5 hours is not a long enough time to see significant change in varying volume and feed. Based on this information, the feed and volume are varied to see the effect on the acetate concentration graph over time. Figures 4 through 6 vary volume while Figures 7 to 9 vary feed. Lastly, Figure 10 starts with ten times the normal initial glucose concentration.
Figure 4. Acetate concentration (g/L) over time for 50 hours. Volume has been changed to be 1.0 L.
Figure 5. Acetate concentration (g/L) over time for 50 hours. Volume has been changed to be 2.5 L.
Figure 6. Acetate concentration (g/L) over time for 50 hours. Volume has been changed to be 3.0 L.
Figure 7. Acetate concentration (g/L) over time for 50 hours. Feed has been changed to be 0.15 g/L.
Figure 8. Acetate concentration (g/L) over time for 50 hours. Feed has been changed to be 0.3 g/L.
Figure 9. Acetate concentration (g/L) over time for 50 hours. Feed has been changed to be 0.5 g/L.
Figure 10. Acetate concentration (g/L) over time for 50 hours. Feed has been changed to be 1.0 g/L.
Figure 11. Acetate concentration (g/L) over time for 50 hours. Initial glucose concentration is multiplied by 10.
As shown in Figures 4 to 6, if the volume is too small, acetate concentration will drop. This is mainly due to the fact that acetate production is related to biomass growth. If the space in which the bacteria grows is too little, the biomass concentration will plateau and acetate production will cease. Acetate concentration also decreases due to the bacteria slowly uptaking acetate as a carbon source. As for the feed, higher values lead to a decrease in acetate concentration over time. Numerically, this is because the negative value in the rate of change is much greater than the production value. This validates that the chosen values for feed, value, and glucose are acceptable.
From the simulation using the feed and volume values displayed in Table 1, the maximum concentration of acetate will be 0.102 g/L. Since the graph starts to plateau around 50 hours, this tells the team that the glucose stock needs to be replenished around that time stamp. This also gives a production rate of 0.00204 g/L/h. While the system could be run for longer, acetate production slows down as time goes on. For an optimized system and a relatively constant acetate production, the system needs to be refreshed before 50 hours. It may also be beneficial to start with a higher initial concentration of biomass, as the acetate concentration at the beginning of the graph is relatively stagnant.
Impact on Wet Lab
The concentration of acetate produced over time gives the wet lab an estimate of the following:
- Ideal volume and feed (V = 2.0L, F = 0.09 g/L)
- Amount of glucose needed
- Acetate production over time
- How often glucose should be injected into the system
- Potential growth-inhibiting stress
Since volume and feed are represented together as a ratio, the ideal volume and feed will be dependent on each other. As found in the analysis by varying both variables, the best values for volume and feed given a 50 hour time frame would be 2.0L for volume and 0.09g/L for feed. Given as a ratio, F/V would be 0.045g. If time were varied, the ratio may have other ideal values.
The initial glucose concentration can be adjusted to generate new acetate concentration graphs over various periods of time. By varying the initial concentration, it can inform the wet lab what concentration they should start with regarding testing. Acetate production rates can be calculated by roughly dividing the total concentration at the end of the period of time by the amount of time the simulation was run for. The acetate concentration over time graphs can also inform when there needs to be more glucose added. When acetate concentration starts to plateau or drop off, that indicates that there is not enough glucose left for the acetate overflow route to take effect. Once the cells use up the glucose, they will begin to use acetate as a carbon source. By looking at the time the acetate concentration starts to plateau, it indicates that more glucose would need to be added.
A more niche but important benefit is upon comparison between the ideal acetate production and the experimental acetate production. Since the acetate production in the simulation is based on Monod-type kinetics, it assumes no stress and perfect growth. There are many sources of stress that could cause the bacteria to produce less than expected. However, the most important factor would probably be if the acetate overproduction puts stress on bacterial growth. If all other sources of stress are eliminated in wet lab and acetate production is still less than normal, the metabolic stress can be quantified based on the difference between ideal and real. This quantification can then be implemented in the simulation to make expectations more accurate.
Future Experiments
Given the scope of the project, time was unfortunately very limited. If we had more time, we would have liked to do some of the following experiments.
It would be a natural experiment to compare the acetate concentration growth to the calculated acetate concentration graph. Since Monod-type kinetics imply ideal growth, stalls in growth and acetate production due to stress would not be predicted. After eliminating as many possible sources of stress as possible, it would have been nice to see if there was some sort of metabolic stress due to acetate overproduction. By averaging the difference in acetate concentrations between ideal and metabolic stress conditions, a mathematical term can be developed and added in to make the model more accurate.
The lethality of acetate concentration was also not explored in this model. There was no mathematical term limiting how high the acetate concentration could go before it would be too concentrated for bacterial growth to occur. With some experimental testing, the concentration of acetate that would be lethal to bacteria can be found. That data could then be brought back into the model to make for more accurate predictions.
An experiment to run would be to check acetate production depending on enzyme concentration. The project utilizes promoters to overexpress the enzymes related to the acetate overflow path. We could then measure the enzyme production levels in the cell. This can then be used to adjust the yield values to get a more accurate model. The reasoning behind going through the extra step of checking based on the enzyme concentration instead of comparing the acetate production directly would be to avoid any stress effects. Getting a theoretical production rate based on the enzyme level would allow the experiment discussed above to be more accurate in predicting metabolic stress.
Summary
For the second module of the wet lab portion of the project, a couple of differential equations are developed and numerically solved to give some insight into the overall system. The ideal volume and feed values are given as 2.0L and 0.09 g/L respectively. The maximum acetate concentration before it starts to plateau is 0.102 g/L and the production rate is given to be 0.00204 g/L/h. Several possible future wet lab experiments are also discussed, which would help make the model more accurate in predicting the acetate concentration over time. These include comparison between ideal growth (from the simulation) and experimental growth, acetate concentration lethality, and enzyme concentration to acetate production correlation.
- Anane, E., López C, D. C., Neubauer, P., & Cruz Bournazou, M. N. (2017). Modelling overflow metabolism in Escherichia coli by acetate cycling. Biochemical Engineering Journal, 125, 23–30. https://doi.org/10.1016/j.bej.2017.05.013
- Averaged Propulsive Body Acceleration (APBA) Can Be Calculated from Biologging Tags That Incorporate Gyroscopes and Accelerometers to Estimate Swimming Speed, Hydrodynamic Drag and Energy Expenditure for Steller Sea Lions. (n.d.). ResearchGate. Retrieved May 27, 2025, from Click Here to View
- Guardia, M. J., & Calvo, E. G. (2001). Modeling of Escherichia coli growth and acetate formation under different operational conditions. Enzyme and Microbial Technology, 29(6), 449–455. https://doi.org/10.1016/S0141-0229(01)00413-6
PHBV Bioreactor
Introduction
While metals are necessary for structural and load-bearing builds, plastic is also an incredibly powerful material for production uses. In aiming to build a colony or station on Mars, it is necessary to utilize resources that already exist there. Plastic, long chains of hydrocarbons, are inherently organic molecules and require carbon sources. Metals can be extracted from the Martian regolith, but the only naturally occurring carbon source on Mars is the atmosphere. The Martian atmosphere is 95% carbon dioxide and an excellent environment for autotrophy.
Our project aims to combine an autotrophic acetogen (acetate producer) and PHB producer in a bioreactor. The autotrophic acetogen would fix carbon from the atmosphere to produce acetate. Acetate can then be fed to the PHB producer for PHB production. Producing a graph of PHB concentration over time not only shows the limit of PHB concentration, but it also allows a prediction of PHB production.
Methods
A system of coupled differential equations describes each organism in the bioreactor. Coupled equations describing biomass, acetate, hydrogen, carbon dioxide concentrations, gas mole balances, and the Monod growth rate are part of the system for an autotrophic acetogen. For the PHB producer, a similar Monod growth rate equation is applied, along with equations measuring rate of change of acetate, nitrogen, and PHB (Cestellos-Blanco, 2021). They are solved numerically through Python code. By recalculating the concentrations and rates of change after each incremental increase in time, graphs of the concentration of acetate and PHB can be produced.
Autotrophic Acetogen Equations
The autotrophic acetogen is similar to the acetate production simulation. The key difference is autotrophy, where the bacteria fixes carbon from atmospheric carbon dioxide. While in individual modules, the acetogen relies on the acetate overflow path and excess glucose. This simulation is based on the WL pathway which combines autotrophy and acetogenic properties in one bacterium (Cestellos-Blanco, 2021). Since the rates of change only rely on the concentrations of the starting material, product, and biomass, this is a valid simulation. Especially once the yield constants are adjusted to fit E. coli and our project’s gene pathway, this would be a good approximation for ideal acetate production.
Biomass (x) is defined using Monod-type kinetics. The change in biomass would be some rate, μ, associated with the growth of bacteria multiplied by the concentration of biomass, denoted using cx. This is shown in Equation 1, below.
Equation 1
The change in concentration of acetate is defined similarly in Equation 2, except with the molar ratio of acetate formed to cells produced (α). This is assuming that the change in acetate concentration will be entirely dependent on biomass growth.
Equation 2
Since the concentration of molecular hydrogen is also consumed in metabolism, it is also an important concentration to take into account. As shown in Equation 3, concentration is dependent on biomass and hydrogen concentrations. The first term in the parentheses indicates the consumption of molecular hydrogen dependent on biomass concentration and bacterial growth rate. It is dependent on yield constants (Y) that define the acetate (Ac) or biomass (X) yield on molecular hydrogen (H2). All of the variables and constants are defined in Table 1. The second term defines how fast molecular hydrogen can dissolve into media to be accessible to bacteria. The difference between the hydrogen saturation concentration and the actual hydrogen concentration is multiplied by kLa, a constant that dictates gas-liquid mass transfer rate.
Equation 3
Equation 4 is very similar to Equation 3 but describing the concentration of carbon dioxide in solution instead of hydrogen.
Equation 4
Equations 5 and 6 control the partial pressure of gaseous hydrogen and carbon dioxide. It is essentially the second term of Equations 3 and 4, describing the mass transfer between gas and liquid phase. It is negated to show that the mass is leaving the gas phase and is thus subtracted from the partial pressure. It is also multiplied by the ratio of the gas and liquid volumes, ideal gas constant, and temperature. The partial pressure needs to be taken into account for a complete closed system because the concentration of the gases will be depleted as the system equilibrates the concentration of gases in media with the partial pressure (p).
Equation 5
Equation 6
The saturation concentration of the gases is dependent on the partial pressure of the atmosphere above the liquid region of the bioreactor. Equations 7 and 8 depict this relation using the Henry constant (H) as a conversion factor.
Equation 7
Equation 8
Lastly, the Monod-type kinetics growth rate is shown in Equation 9, as a function of concentration of hydrogen and carbon dioxide.
Equation 9
PHB Producer Equations
The change in biomass concentration (Equation 10) is modeled similarly to the one for autotrophic acetogen (Equation 1), except that the growth rate (μx) is defined using acetate and nitrogen.
Equation 10
Similarly to Equation 2, the production rate of PHB (Equation 11) is associated with the biomass concentration. This is also assuming that the production of PHB is dependent on bacterial growth.
Equation 11
Acetate is the precursor to PHB and is thus consumed as the PHB producer concentration grows and PHB is produced, as shown in Equation 12. The first term describes acetate concentration decreasing due to acetate being the carbon source of bacterial growth. The second term refers to the consumption of acetate due to PHB production.
Equation 12
Since PHB production doesn’t depend on nitrogen, Equation 13 only includes the first term from Equation 12, adjusted for nitrogen instead of acetate.
Equation 13
The rate at which biomass increases (Equation 14) is similar to the Monod-type kinetics growth rate in Equation 9. However, the variables are changed to depend on acetate and nitrogen instead, and include an inhibition term.
Equation 14
Lastly, the rate of PHB production is shown in Equation 15. The variable fPHB is the PHB to biomass ratio. It is corrected by raising it to the power of some β. This term limits the PHB concentration so that the concentration of PHB does not reach lethal levels.
Equation 15
Variables and Definitions
The variables for each of the equations are defined in the following tables, along with the associated value for E. coli if applicable. The bolded number in the Equation column is the equation that the variable is defined in if it is a function. The tables are organized by chronological appearance in the equations. Table 1 are all the variables for the autotrophic acetogens. Table 2 are all the variables for the PHB producers.
Table 1. List of variables organized by equation for the autotrophic acetogen (E.coli). The bolded number is the equation that the variable is defined in. *See Assumptions section.
Table 2. List of variables organized by equation for the PHB producer (E. coli). The bolded number is the equation that the variable is defined in. *See Assumptions section.
These equations can then be numerically calculated in Python. The simulation starts at time equals 0, and then adds on small steps of time. After each small step of time is added, the variables that are listed as “continuously calculated” are then recalculated. Since the rate of change is then known, it can be used to calculate the new concentrations. Using these concentrations, the concentration for biomass, acetate, hydrogen, CO2, and PHB can be graphed over time.
Assumptions
As this simulation is a model of the real system, a couple assumptions are made for ease of calculation:
- Monod-type kinetics growth
- Excess dissolved oxygen (DOT)
- Some K constants would be approximately the same for E. coli as in the Cestellos-Blanco (2021) paper
The growth curves are written with Monod-type kinetics in mind. This assumes that there is no metabolic stress due to acetate or PHB production, or any unexpected stress that may come with being in a bioreactor.
A couple of the values are educated guesses. In adjusting for element analysis, a stoichiometric model of E. coli was referenced, resulting in 25.62 gCDW/mol (Pramanik 1997). Since this is within significant figures of the 25 gCDW/mol referenced in Cestellos-Blanco, 2021, the element balance values are taken to be accurate enough for the purposes of this model. This simulation has also assumed excess dissolved oxygen, assuming that the bioreactor would be somewhere where the partial pressure of oxygen is enough for a long runtime.
Another assumption is that the half-velocity constants in Table 1 and the constants associated with nitrogen in Table 2 are approximately the same for E. coli and the bacterium used in Cestellos-Blanco, 2021. Another assumption is that the difference in values will not have a significant effect on the shape of the curve or production time.
Results
Before making the curves using the E. coli values, the simulation is tested using the values from Cestellos-Blanco, 2021. This is to ensure that the initial simulation results are reasonable and that there were no problems with the code and underlying differential equations. The simulation was run for 27 and 28 hours to identify when the system would fail. As seen in Figure 1 (27 hours) and Figure 2 (28 hours), PHB production peaks at around 25 hours and then equilibrates until the system fails at 28 hours. This is because the bacteria will begin to use PHB as a carbon source in the event that biomass concentration becomes too large. The system failing means that the concentration for PHB immediately falls off to be 0, as seen in Figure 2.
Figure 1. PHB concentration vs time using the values from Cestellos-Blanco, 2021, run for 27 hours.
Figure 2. PHB concentration vs time using the values from Cestellos-Blanco, 2021, run for 28 hours.
Knowing that the system does fail eventually, the simulation with E. coli values was run until failure (Figure 3 and 4). This simulation fails at 41 hours, as seen in Figure 4. It is interesting to note that while the E. coli system does equilibrate, it peaks at a PHB concentration slightly higher than the equilibrium. The maximum concentration is 158 mM at 36 hours. The production rate is also roughly exponential, which is natural given that biomass growth is also exponential and that PHB concentration rate of change is dependent on biomass growth.
Figure 3. PHB concentration vs time using E. coli values, run for 40 hours.
Impact on Policy and Practice
Since production of PHB is a priority in this bioreactor, it would be best if the bioreactor were refreshed at the peak of PHB concentration. That way, the maximum amount of PHB would be extracted at a given instance. If this bioreactor were to be implemented, PHB extraction happens every 36 hours, for a concentration of 158mM. In other words, every three days can produce 316mM of PHB. This provides Policy and Practice with two important pieces of information:
- PHB production rate
- Labor cost
In calculating the cost of a space mission, all time and resources are accounted for. Having as much data as possible for an accurate estimate is critical for cost efficiency. Plastic is a powerful material for producing everyday objects. The ease of plastic production and versatility of applications is what makes it such a widespread material today. While it is hugely detrimental to the environment, the advantages of plastic material cannot be understated. Luckily, the added benefit of PHB being a bioplastic makes it easy to degrade and reuse carbon components in a closed loop system.
Accounting for PHB production rate over large time scales gives reliable timelines for establishing a functioning colony. For example, some pieces of mission critical equipment may require a certain amount of plastic. Knowing production rate, the timeline for producing that amount of material can be estimated. Since the production rate is given in concentration per hour, the volume of the bioreactor can be scaled up or down for efficiency.
The labor cost of maintaining a bioreactor can also be accounted for. Astronaut time is one of the more difficult things to budget, as it is not solid and quantifiable. On a mission, there are so many experiments and data collection to do that it is very necessary to budget astronaut working time. Knowing that the bioreactor would have to be maintained every 36 hours means that an astronaut would have to dedicate some time roughly every other day to refresh the bioreactor and extract PHB, unless some mechanical component would be able to do it automatically.
Future Experiments
As an idealized model of the project, there are a couple experiments that could be run to improve the accuracy of the simulation.
The main experiment would be to produce Michaelis-Menten kinetic graphs in order to get accurate values for the assumed K variables. This would require measuring the rate of reaction as the substrate concentration (hydrogen, CO2, and nitrogen) increases. KS,H2, KS,CO2, and KS,N can then be calculated by finding substrate concentration when the velocity (reaction rate) is half of the maximum velocity. Not only that, several experiments validating the element analysis-derived yield values should also be run.
The concentration of PHB over time should also be taken and compared to the simulation. As with the acetate simulation, Monod-type kinetics growth assumes no metabolic stress or environmental stress-inducers. If there is a difference in PHB concentration, the stress can be quantified and introduced to the simulation for more accurate estimates.
Summary
The idealized version of the PHB bioreactor is modeled for the evolution of PHB concentration over time. Differential equations are developed to describe the rate of change of biomass, acetate, PHB, and starting material over time. These equations are solved numerically using Python. The simulation can be run for 40 hours before it fails, and peaks at a concentration of 158mM at 36 hours. These values are useful for calculating mission costs such as scaling the size of bioreactor or labor costs. Some future wet lab experiments for improving the model are discussed, for validation of the value of some constants used in simulation.
- Andersen, K. B., & von Meyenburg, K. (1980). Are growth rates of Escherichia coli in batch cultures limited by respiration? Journal of Bacteriology, 144(1), 114–123. https://doi.org/10.1128/jb.144.1.114-123.1980
- Cestellos-Blanco, S., Friedline, S., Sander, K. B., Abel, A. J., Kim, J. M., Clark, D. S., Arkin, A. P., & Yang, P. (2021). Production of PHB From CO2-Derived Acetate With Minimal Processing Assessed for Space Biomanufacturing. Frontiers in Microbiology, 12. https://doi.org/10.3389/fmicb.2021.700010
- CRC Handbook of Chemistry and Physics, 84th edition. (CRC Press, 2004). doi:10.1136/oem.53.7.504
- de Oliveira Santos, W., David de Oliveira, R., Cabrera Gomez, J. G., & Antonio Carrillo Le Roux, G. (2025). Analysis of Poly-3-Hydroxybutyrate Production with Different Microorganisms Using the Dynamic Simulations for Evaluation of Economic Potential Approach. ACS Omega, 10(26), 27756–27774. https://doi.org/10.1021/acsomega.4c11178
- Dreybrodt, W., Lauckner, J., Zaihua, L., Svensson, U., & Buhmann, D. (1996). The kinetics of the reaction CO2 + H2O → H+ + HCO3− as one of the rate limiting steps for the dissolution of calcite in the system H2OCO2CaCO3. Geochimica et Cosmochimica Acta, 60(18), 3375–3381. https://doi.org/10.1016/0016-7037(96)00181-0
- Dry weight of E. coli B/r cell in exponential—Bacteria Escherichia coli—BNID 103904. (n.d.). Retrieved August 27, 2025, from https://bionumbers.hms.harvard.edu/bionumber.aspx?id=103904&ver=17
- Dukes, A. D. (2020). Measuring the Henry’s Law Constant for Carbon Dioxide and Water with UV-visible Absorption Spectroscopy. Analytical Sciences: The International Journal of the Japan Society for Analytical Chemistry, 36(8), 971–975. https://doi.org/10.2116/analsci.19P477
- Haberle, R. M. (2015). SOLAR SYSTEM/SUN, ATMOSPHERES, EVOLUTION OF ATMOSPHERES | Planetary Atmospheres: Mars. In G. R. North, J. Pyle, & F. Zhang (Eds.), Encyclopedia of Atmospheric Sciences (Second Edition) (pp. 168–177). Academic Press. https://doi.org/10.1016/B978-0-12-382225-3.00312-1
- Liu, Z., Lin, B., Campbell, J. F., Yu, J., Geng, J., & Jiang, S. (2024). Martian column CO2 and pressure measurement with spaceborne differential absorption lidar at 1.96 µm. Atmospheric Measurement Techniques, 17(9), 2977–2990. https://doi.org/10.5194/amt-17-2977-2024
- Lucile, F., Cézac, P., Contamine, F., Serin, J.-P., Houssin, D., & Arpentinier, P. (2012). Solubility of Carbon Dioxide in Water and Aqueous Solution Containing Sodium Hydroxide at Temperatures from (293.15 to 393.15) K and Pressure up to 5 MPa: Experimental Measurements. Journal of Chemical & Engineering Data, 57(3), 784–789. https://doi.org/10.5194/amt-17-2977-2024
- Miu, D.-M., Eremia, M. C., & Moscovici, M. (2022). Polyhydroxyalkanoates (PHAs) as Biomaterials in Tissue Engineering: Production, Isolation, Characterization. Materials, 15(4), 1410. https://doi.org/10.3390/ma15041410
- Mozumder, S., Garcia-Gonzalez, L., De Wever, H., & Volcke, E. (2015). Poly(3-hydroxybutyrate) (PHB) production from CO2: Model development and process optimization. Biochemical Engineering Journal, 98. https://doi.org/10.1016/j.bej.2015.02.031
- PHB granules. – Biodegradable for Sustainable Packaging. (n.d.). https://www.goodfellow.com/usa/phb-granules-group#properties
- Pramanik, J., & Keasling, J. D. (1997). Stoichiometric model of Escherichia coli metabolism: Incorporation of growth-rate dependent biomass composition and mechanistic energy requirements. Biotechnology and Bioengineering, 56(4), 398–421. https://doi.org/10.1002/(SICI)1097-0290(19971120)56:4<398::AID-BIT6>3.0.CO;2-J
- Xiao, Y., Feng, X., Varman, A. M., He, L., Yu, H., & Tang, Y. J. (2012). Kinetic Modeling and Isotopic Investigation of Isobutanol Fermentation by Two Engineered Escherichia coli Strains. Industrial & Engineering Chemistry Research, 51(49), 15855–15863. https://doi.org/10.1021/ie202936t
PHBV Properties
Introduction
While PHB itself, a homopolymer, is relatively rigid and brittle, PHBV, a copolymer of 3-hydroxybutyrate (3HB) and 3-hydroxyvalerate (3HV) monomers, is a stronger and more flexible composite material. Knowing the properties of different PHBV compositions, based upon their different HB:HV ratios, can help make predictions about which constitutions of PHBV to use for different manufacturing applications.
The positive association between 3HV composition and both durability and flexibility of material has been described in material science literature. One article in particular used the parameters of tensile strength, elongation, and elastic modulus to characterize the physical properties of materials made with different percentages of 3HV.
Tensile strength, measured here in megapascals, can be defined as the maximum amount of stress a certain material can withstand while being stretched before breaking. Elongation describes how much the material can be stretched before this breaking point, in this case expressed as a percentage of the material’s original length. Elastic modulus, measured here in gigapascals, defines the stiffness or resistance to deformation of a material being stretched. It is calculated by taking the ratio of tensile stress (applied force per unit area) to strain (resulting deformation). Thus, a higher elastic modulus will indicate a more rigid material.
Methods
A linear regression approximation can be made based on the experimental values. Graphing the correlation between 3HV percentage and tensile strength, elongation, and elastic modulus produces an equation for each of the material properties. The experimental data is inputted into MS Excel, graphed, and fitted to linear regression models.
Results
Graphing the empirical data from Avella, 2000, a best fit line can be generated, showing a roughly linear relationship between material properties and 3HV percent composition (Figures 1, 2, and 3).
Figure 1. Linear regression of tensile strength of PHBV dependent on 3HV percent composition.
Figure 2. Linear regression of elongation of PHBV dependent on 3HV percent composition.
Figure 3. Linear regression of elastic modulus of PHBV dependent on 3HV percent composition.
As shown in Figures 1, 2, and 3, as the composition of HV within the PHBV copolymer increases, the tensile strength of the PHBV decreases, its elongation percentage increases, and its elastic modulus decreases.This can be interpreted to mean that the greater the percentage 3HV within a material, the less resistant it is to deformation resulting from stress and thus the more flexible it will be under pressure. This relationship can be taken into account when deciding what materials will need to withstand or flex under pressure. Equations for each linear regression are shown in Figures 1, 2, and 3, and also listed in Table 1.
Table 1. R2 values and equations of each linear regression model.
As seen in Table 1, the R2 values indicate that a linear regression model is appropriate for the data set. The absolute values of all three of the R2 values are above 0.8, validating the choice of using a linear regression model.
Implications
The linear regression equations produced are incredibly helpful for developing PHBV with a certain elasticity or tensile strength end goal. In terms of our project, changing the PHBV would look like tweaking the threonine concentrations. In wet lab module 4, tdcB and tdcE (threonine dehydratase and 2-ketobutyrate formate-lyase) convert threonine into propionyl-CoA, a precursor to 3HV. Using the equations to find that there needs to be some mol% of 3HV, more or less threonine can be added. More threonine would increase the concentration of 3HV and consequently the mol% of 3HV in the final PHBV product.
Future Experiments
An experiment to relate our project’s wet lab component to PHBV percent composition is to test for the minimum and maximum percent composition. Since threonine is naturally made, there will always be some baseline 3HV in the product. The percent composition of 3HV when there is no added threonine would be the minimum. It is infeasible to produce anything less than this lower limit because it is impossible to remove naturally occurring threonine. For the maximum percent composition, this can be found empirically by checking the PHBV percent compositions of the bacteria grown in varying concentrations of excess threonine. Injecting threonine should increase the 3HV content in PHBV, since more of the 3HC precursor is being made. However, there should be a limit as to how high the mol% of 3HV can be. 3HB is consistently being made regardless of threonine concentration. This means that when the enzymes producing propionyl-CoA are at maximum velocity, there will still be 3HB in solution. At maximum velocity, the PHBV produced would have the maximum mol% 3HV. Experimentally, this would mean that despite threonine concentration increasing, percent composition would plateau. Knowing the minimum and maximum mol% 3HV would allow a set range of 3HV percent compositions to be determined.
Summary
After fitting the data to a linear regression, equations were developed and the R2 values were obtained. Based on the R2 values, a linear model is appropriate for the data and linear equations can be reliably used to predict material properties for given mol% 3HV. The reverse can also be true, as having an end goal tensile strength, elongation, or elastic modulus can yield an approximate %mol 3HV in the final PHBV product. Future experiments are also discussed to find the limits on the range of %mol 3HV.
- Avella, M., Martuscelli, E., & Raimo, M. (2000). Review Properties of blends and composites based on poly(3-hydroxy)butyrate (PHB) and poly(3-hydroxybutyrate-hydroxyvalerate) (PHBV) copolymers. Journal of Materials Science, 35(3), 523–545. https://doi.org/10.1023/A:1004740522751
- Mazur, K. E., Jakubowska, P., Gaweł, A., & Kuciel, S. (2022). Mechanical, thermal and hydrodegradation behavior of poly (3-hydroxybutyrate-co-3-hydroxyvalerate) (PHBV) composites with agricultural fibers as reinforcing fillers. Sustainable Materials and Technologies, 31, e00390. https://doi.org/10.1016/j.susmat.2022.e00390
