Figure 1 | Mathematical Modeling Framework

In the mathematical modeling section, we explored the role of the cascade system and established a microbial interaction analysis system based on growth kinetics and statistical hypothesis testing. By comparing the biomass, area under the curve, and maximum specific growth rate under mono and co-culture conditions, we identified competitive, neutral, or reciprocal relationships among strains and integrated them into a system for the convenience of other researchers. In addition, the pigment diffusion and adsorption processes were described by using partial differential equations, achieving quantitative prediction and system optimization of staining behavior.

Due to the difficulty in quantitatively observing some processes of biological systems and the difficulty in accurately predicting their dynamic responses, the verification of experimental feasibility is of vital importance. In this study, an ordinary differential equation model based on biological expression was constructed using Python to simulate and analyze biological processes to explore the effect of cascade amplification and verify the design of experimental pathways. This not only provides scientific guidance for wet experiments but also ensures the reliability of the design, reduces research costs and time.

Table 1 | Symbol explanation

Symbol	Explanation
$M_\text{ECF}$	The mRNA concentration of ECF
$P_\text{ECF}$	ECF protein concentration
$M_\text{LuxI}$	The mRNA concentration of LuxI
$\delta_m$	mRNA degradation rate
$\delta_p$	Protein degradation rate
$K_1$	The Michaelis constant of transcription
$K_2$	The Michaelis constant for translation

To explore the cascade amplification effect of the downstream luxI protein on ECF expression and cellulose synthesis in the pathway, we constructed a transcription-translation dynamics model based on ordinary differential equations and compared and analyzed the differences in the system with and without the downstream luxI protein. In this model, both the transcription and translation processes follow the Michaelis-Menten equation, that is, the reaction rate is in a saturation relationship with the substrate concentration.

The model assumes that the upstream components can be continuously expressed, and the downstream luxI protein can further activate the promoter to enhance the transcription rate of the ECF protein, thereby accelerating the subsequent generation of cellulose. The corresponding differential equation of dynamics is as follows: $$\begin{cases}\displaystyle\frac{dP_\text{LuxI}}{dt}=\frac{V_1M_\text{LuxI}}{K_1+M_\text{LuxI}}-\delta_pP_\text{LuxI}\\[3ex] \displaystyle\frac{dM_\text{ECF}}{dt}=\frac{V_2P_\text{LuxI}}{K_2+P_\text{LuxI}}-\delta_mM_\text{ECF}\\[3ex] \displaystyle\frac{dP_\text{ECF}}{dt}=\frac{V_1M_\text{ECF}}{K_1+M_\text{ECF}}-\delta_pP_\text{ECF}\\ \end{cases}$$

Based on the "transcription-translation-degradation" kinetic law of gene expression, we wrote simulation code in Python to quantitatively compare the dynamic characteristics of ECF protein concentration under different conditions. The comparison result generated by the code running is shown in the figure:

Figure 2 | Compare the ECF Production with and without Downstream luxI Protein

The results show that in the design system containing downstream luxI protein, the accumulation of ECF protein within the same period of time is significantly higher than that in the design without luxI protein. This difference provides favorable conditions for the subsequent efficient generation of cellulose.

In the actual experimental process, factors such as light and temperature can have an impact on the reaction. To prove that our reaction can proceed stably under typical and different environmental conditions, we introduced disturbance factors in the simulation process. These disturbances range from -15% to 15%, with seven different levels. Through our simulation, the trend of ECF protein changes was obtained, as shown in the figure:

Figure 3 | Concentrations of Substances in Different Disturbed Systems

After analysis, it was found that under different disturbance conditions, the generation trends of ECF proteins were similar with small errors. This indicates that our reaction has good stability and can operate effectively in different environments.

Table 2 | Symbol Description

Symbol	Explanation
$K$	Environmental carrying capacity
$AUC$	Area under the growth curve
$\mu_\text{max}$	Maximum specific growth rate
$II$	Interaction Index
$\alpha,\beta$	Competition coefficients
$R$	Initial inoculation ratio of Strain A to B
$t_s$	Temperature switch time
$Y_P,Y_Q$	Cellulose yield coefficient, Pigment yield coefficient
$P,Q$	Cellulose yield, Pigment yield
$X_A,X_B$	Biomass of Strain A, Biomass of Strain B

Analysis of Interaction Types in Microbial Co-Culture

Our project involves two types of E.coli : Strain A (cellulose producer) and Strain B (pigment producer). When they are cultivated in the same system, they may compete for resources, or they may not interfere with each other, or be mutually beneficial under specific conditions. Therefore, we aim to analyze the types of interactions between the two by comparing the growth kinetics parameters of their monoculture and co-culture.

Collect Data and Calculate Growth Parameters

To accurately assess the interaction, we communicated with the experimental group and set up three culture conditions to obtain their growth curves:

• Single-cell culture: Strain A is cultured alone in the culture medium.

• Single-cell culture: Strain B is cultured alone in the culture medium.

• Co-culture group: Co-culture strain A and B in an initial ratio of 1:1.

For each growth curve, three key parameters are extracted:

• Maximum biomass K: The carrying capacity of a population in a specific environment;

• Area under the curve (AUC) : To measure the overall growth performance, the trapezoidal rule is adopted for calculation: $$\text{AUC}=\sum_{i=0}^{n-1}\frac{\text{OD}_i+\text{OD}_{i+1}}2·(t_{i+1}-t_i)$$

• Maximum specific growth rate $\mu_\text{max}$ : The fastest proliferation capacity of microorganisms under ideal conditions: $$\ln(\text{OD})=\mu_\text{max}·t+b$$ Where, $\mu_\text{max}$ ($h^{-1}$) represents the slope of the fitted straight line, and its value directly characterizes the inherent growth potential of the strain in a specific environment.

Statistical Hypothesis Testing

The interaction types between the two are determined respectively from the total effect test and the individual effect test, and the interaction index is introduced.

1.Total Effect Test

$H_0$ : There is no significant difference between co-culture and individual culture, that is, the interaction is "neutral".

$H_1$ : There are significant differences between co-culture and individual culture, that is, there exists a "competitive" or "reciprocal" relationship.

Calculate the predicted value, representing the theoretical expectation under the "neutral" relationship:

• Predicted maximum biomass $K_\text{predicted}$ : $$K_\text{predicted}=K^\text{solo}_A+K^\text{solo}_B$$

• Predict the area under the curve $\text{AUC}_\text{predicted}$ : $$\text{AUC}_\text{predicted}=\text{AUC}^\text{solo}_A+\text{AUC}^\text{solo}_B$$

Finally, conduct a statistical test. The one-sample t-test is used to compare whether there is a significant difference between the mean and the predicted value of a set of data. The calculation formula is: $$t=\frac{\overline x-\mu}{s/\sqrt n}$$ P value < 0.05 indicates that the null hypothesis is rejected at the significance level of 0.05($\alpha=0.05$).

2.Individual Effect Test

To avoid masking asymmetric interactions by relying only on total biomass or AUC, we introduced an individual effect test, analyzing co-culture impacts per strain with paired t-tests on mono- and co-culture data to reduce individual variation and enhance statistical sensitivity. $H_0$ and $H_1$ are the same as above, and the statistical test quantity is: $$t=\frac{\overline d-\mu_0}{s_d/\sqrt n}$$ Where, $d$ is the paired difference value: $$d_{A_i}=K^\text{co}_{A_i}-K^\text{solo}_{A_i},\quad d_{B_i}=K^\text{co}_{B_i}-K^\text{solo}_{B_i}$$

3.Interaction Index

To make up for the fact that significance tests can only determine "whether" but not assess "the degree", an interaction index is introduced.

Interaction index based on maximum biomass and the area under the curve: $$II_K=\frac{K_\text{co}-K_\text{predicted}}{K_\text{predicted}}$$ Interaction index based on the area under the curve: $$II_\text{AUC}=\frac{\text{AUC}_\text{co}-\text{AUC}_\text{predicted}}{\text{AUC}_\text{predicted}}$$

Criteria for Judging the Type of Interaction

Table 3 | Criteria for Judging Interaction Types

Total effect test	Interaction index	Individual effect test	Type
$p < 0.05$	$II < -0.2$	The biomass of both strains decreased significantly.	Strong competition
	$-0.2 \leq II < -0.1$	At least one strain decreased significantly, while the other strain remained unchanged or decreased.	Moderate competition
	$-0.1 \leq II < -0.05$	One strain decreased significantly, while the other showed no significant change.	Weak competition/bias
	$II > 0.2$	The biomass of both strains increased significantly.	Strong reciprocity
	$0.05 < II \leq 0.1$	At least one strain increased significantly, while the other strain remained unchanged or increased.	Moderate reciprocity
$p \geq 0.05$	$-0.05 \leq II \leq 0.05$	Neither of the two strains showed significant changes.	Neutral
	$II < -0.1, II> 0.1$	The individual test results are inconsistent or vary greatly.	Uncertain/Complex

Definition and Calculation of Competitiveness

After determining the type of interaction, we want to further quantify the intensity of the interaction to facilitate subsequent modeling.

Situation 1: Competitive Relationship

Based on the assumption that "the population growth rate is zero in a stable competitive state", the Lotka-Volterra model can derive the competition coefficient that quantifies the intensity of mutual inhibition among species.

For strain A and strain B , their competition coefficients are respectively defined as: $$\alpha=\frac{K^\text{solo}_{A}-K^\text{co}_{A}}{K^\text{co}_{B}},\quad\beta=\frac{K^\text{solo}_{B}-K^\text{co}_{B}}{K^\text{co}_{A}}$$

Where, $\alpha$ represents the competition coefficient of strain B over strain A. The larger the $\alpha$ value is, the stronger the inhibitory ability of strain B on strain A is. $\beta$ is the same as above. $K_A^\text{solo},K_B^\text{solo}$ are the environmental carrying capacity of strains A and B. $K_A^\text{co},K_B^\text{co}$ is the biomass of strains A and B when they reach a stable state through co-culture.

Thus, the competitiveness index CI is obtained: $$CI=\frac\beta\alpha$$

Situation 2: Neutral Relationship

In a neutral relationship, the strain that grows faster can utilize resources earlier and gain an advantage in the early stage of population establishment. Define the relative growth competitiveness RC: $$RC=\frac{r_{A,\max}}{r_{B,\max}}$$

Where, $r_{A,\text{max}},r_{B,\text{max}}$ represent the maximum specific growth rate obtained by A and B under separate culture.

Situation 3: Reciprocal Relationship

Cooperative efficiency is defined as the ratio of the total output of the co-culture system to that of the individual culture system: $$CE=\frac{K_{\text{total}}^{\text{co}}}{K_{\text{A}}^{\text{solo}}+K_{\text{B}}^{\text{solo}}}$$

Where, $K_{\text{total}}^{\text{co}}=K_{\text{A}}^{\text{co}}+K_{\text{B}}^{\text{co}}$ represents the total biomass of co-culture.

Determine the Optimal Proportion of Application and the Time for Temperature Switching

Based on the analysis of strain interactions and competitiveness, a model is established to determine key production parameters: the initial inoculation ratio $R$ and the temperature switching time $t_{\text{switch}}$.

The objective function is the maximum yield of cellulose and pigment: $$\text{Obj}(R,t_{\text{switch}})=w_1·P_\text{cellulose}(t_{\text{switch}})+w_2·P_\text{pigment}(t_{\text{end}})$$ Where, $P_\text{cellulose}(t_{\text{switch}})$ represents the cellulose yield accumulated at time $t_{\text{switch}}$ ; $P_\text{pigment}(t_{\text{end}}$ represents the accumulated pigment output at time $t_{\text{end}}$ at the end of production, and $w_1,w_2$ is the weight factor).

1. Obtain Key Parameters

$Y_{P1/X1}$ represents the cellulose yield coefficient. $Y_{P2/X2}$ represents the pigment yield coefficient.

$r_{A25},r_{B25},r_{B37}$ represents the maximum specific growth rates of strain A at 25°C and Strain B at 25°C and 37°C respectively.

2. Predict the biomass of the bacteria

Assuming the total initial inoculation density is fixed. R is the initial inoculation ratio. $D(0)$ represents the total initial biomass. $$X_A(0)=\frac{R}{R+1}\times X(0),\quad X_B(0)=\frac{1}{R+1}\times X(0)$$

Where, $D(0)$ represents the total initial biomass.

25°C

If there is a competitive relationship between the two strains, the Lotka-Volterra competitive model is used to describe their dynamic process, and the equations are as follows:

$$ \begin{cases}\displaystyle \frac{dX_A}{dt}=r_{A25}·X_A(1-\frac{X_A+\alpha_{AB}·X_B}{K_A})\times f_{\text{comp}}(A,B)\\[3ex] \displaystyle \frac{dX_B}{dt}=r_{B25}·X_B(1-\frac{X_B+\alpha_{BA}·X_A}{K_B})\times f_{\text{comp}}(A,B) \end{cases} $$

As the cultivation process continues, the intensity of interaction between the strains will change. Therefore, a resource competition adjustment function is added:

$$f_{\text{comp}}(A,B)=\frac{1}{1+\exp(-k_{\text{resource}}\times(S_0-\rho_AX_A-\rho_BX_B))}$$

Among them, $S_0$ represents the initial substrate concentration, $\rho_A,\rho_B$ are the resource consumption rate per unit biomass, and $k_{\text{resource}}$ is the resource sensitivity parameter.

If there is a neutral relationship between the two strains, they grow independently without affecting each other:

$$\begin{cases} \displaystyle \frac{dX_A}{dt}=r_{A25}·X_A(1-\frac{X_A}{K_A})\\[3ex] \displaystyle \frac{dX_B}{dt}=r_{B25}·X_B(1-\frac{X_B}{K_B})\\ \end{cases}$$

If there is a reciprocal relationship between the two strains, the existence of the other population promotes the growth of this population: $$\begin{cases} \displaystyle \frac{dX_A}{dt}=r_{A25}·X_A(1-\frac{X_A}{K_A+\beta_{AB}·X_B})\times g_\text{metab}(A,B)\\[3ex] \displaystyle \frac{dX_B}{dt}=r_{B25}·X_B(1-\frac{X_B}{K_B+\beta_{AB}·X_A})\times g_\text{metab}(A,B)\\ \end{cases}$$

The metabolites produced during the culture process can also affect the growth of the strain. Add the metabolite exchange regulation function: $$g_\text{metab}(A,B)=1+\eta[1-\exp(-\lambda·\min\{A,B\})$$ Where, $\eta$ represents the maximum metabolic promotion intensity, and $\lambda$ represents the diffusion efficiency of metabolites.

37°C

The growth of strain A is strongly inhibited, while strain B grows rapidly: $$\begin{cases} \displaystyle \frac{dX_A}{dt}=r_{A37}·X_A·h_\text{inhibit}(t)-\delta_AX_A\\[3ex] \displaystyle \frac{dX_B}{dt}=r_{B37}·X_B(1-\frac{X_B}{K_B})\times [1+\omega\frac{X_A}{K_A}]\\ \end{cases}$$

Where, $h_\text{inhibit}$ is the temperature inhibition function, indicating that bacteria A are gradually inhibited as the temperature rises: $$h_\text{inhibit}=1/[1+\exp(k_\text{temp}(T-T_\text{critical})]$$ Where, $T_\text{critical}$ is the critical growth temperature of strain A, $k_\text{temp}$ is the temperature inhibition intensity, $\delta_A$ is the decay rate of strain A at 37°C, and $\omega$ is the residual effect coefficient of the residual strain A on strain B.

3. Predict Product Output

Cellulose yield: $$P_\text{cellulose}(t_{\text{switch}})=Y_{P1/X1}·(X_A(t_{\text{switch}})-X_A(0))$$ Pigment production: Mainly produced by strain B at 37°C. $$P_\text{pigment}(t_{\text{end}})=Y_{P2/X2}·(X_B(t_{\text{end}})-X_B(t_{\text{switch}})$$

4. Grid Search Algorithm

The optimal result is obtained by using the Grid Search algorithm. The steps are as follows:

1. Define the parameter space: $R\subset[0.1,10]$, $t_{\text{switch}}\in[5,25]$ hours.

2. Generate parameter mesh: Take points uniformly within the parameter range.

3. Traversal and calculation: Calculate the $P_\text{cellulose}$ and $P_\text{pigment}$ of each group of parameters and calculate the value of $\text{Obj}$.

4. Seek the optimal solution: The parameters corresponding to when $\text{Obj}$ reaches its maximum value.

Interactive Applications

To simplify the analysis of microbial interactions, we have developed an interactive desktop application integrating data computing and visualization based on Python PyQt5. This application supports a series of functions from CSV data upload and verification to interaction type determination, degree quantification, parameter optimization and growth curve drawing, providing an efficient and scientific computing platform for the research and industrialization of microbial co-culture.

To verify the accuracy of the system, we imported the data of the experimental group into this interactive system for analysis. As shown in the figure, the growth dynamics and interaction patterns of the two bacterial strains output by the system are highly consistent with the known biological characteristics, thereby verifying the reliability of the system in the study of microbial interactions.

Figure 4 | Interface of Interaction Analysis Results

Figure 5 | Showing the Interface of Parameter Optimization Results

For specific operation and program download methods, please refer to the "User Guide for Microbial Co-culture Interaction Analysis System".

Our one-step biosynthesis, the project's core innovation, occurs within the porous BC matrix, producing pigments and cellulose simultaneously. Unlike conventional ODE-based dyeing models, this process requires a PDE framework to capture spatiotemporal dye transport and adsorption alongside temporal cellulose-to-pigment dynamics.

PDE Model

To model dye migration and adsorption in the bacterial cellulose (BC) membrane, we consider a representative elementary volume (REV) with liquid-filled pores and solid matrix. Solute concentration in pores is $C(x,t)$,and adsorption on the solid is $q(x,t)$. Without external reactions, total solute accumulation balances the divergence of diffusive flux: $$\phi\frac{\partial C}{\partial t}+\frac{\partial q}{\partial t}=\nabla\cdot\bigl(D(C)\nabla C\bigr)$$

Here, $\phi$ is porosity and $D(C)$ the effective diffusion coefficient. By referencing solute concentration to pore volume ($\phi = 1$), the equation simplifies to: $$\frac{\partial C}{\partial t}=\nabla\cdot\bigl(D(C)\nabla C\bigr)-\frac{\partial q}{\partial t}​$$

The diffusive flux follows Fick's law, $\mathbf{J}=-D(C)\nabla C$, with $D(C)$ typically decreasing with concentration. For modeling, an empirical Michaelis–Menten–type expression is used: $$\quad D(C)=\frac{D_0}{1+\alpha C}$$

Here, $D_0​$ denotes the diffusion coefficient in the dilute limit, while $\alpha$ is the concentration-inhibition parameter, reflecting intermolecular interactions or pore-blocking effects.

Adsorption kinetics at the solid–liquid interface follow a Langmuir-type model, with adsorption proportional to solute concentration and available sites, and desorption proportional to occupied sites, yielding: $$\frac{\partial q}{\partial t}=k_a C\,(q_{\max}-q)-k_d q$$

Here, ⁡$q_{\max}$​ is the saturation capacity, and $k_a​$ , $k_d$ are adsorption/desorption rate constants. At equilibrium, this reduces to the Langmuir isotherm, leading to the coupled single-continuum model: $$\begin{cases} \displaystyle \frac{\partial C}{\partial t}=\nabla\cdot\bigl(D(C)\nabla C\bigr)-\frac{\partial q}{\partial t}\\[3ex] \displaystyle \frac{\partial q}{\partial t}=k_a C\,(q_{\max}-q)-k_d q\end{cases}​$$

To model BC's dual-porosity, the system is split into macropores ($\phi_m​$), and and micropores ($\phi_s​$), with $\phi_m + \phi_s + \phi_{\text{solid}} = 1$ . Solute concentrations $C_m$​ and $C_s$​ diffuse with $D_m(C_m)$ , and mass exchange is modeled by $k_{ex}(C_m - C_s)$ , yielding: $$\begin{cases} \displaystyle \phi_m\frac{\partial C_m}{\partial t}=\nabla\cdot\bigl(D_m(C_m)\nabla C_m\bigr)-\phi_m\frac{\partial q_m}{\partial t}-k_{ex}(C_m-C_s)\\[3ex] \displaystyle \phi_s\frac{\partial C_s}{\partial t}=\nabla\cdot\bigl(D_s(C_s)\nabla C_s\bigr)-\phi_s\frac{\partial q_s}{\partial t}+k_{ex}(C_m-C_s) \end{cases}​$$

The corresponding adsorption kinetics are given by $$\frac{\partial q_m}{\partial t}=k_aC_m\,(q_{m,\max}-q_m)-k_dq_m, \qquad \frac{\partial q_s}{\partial t}=k_aC_s\,(q_{s,\max}-q_s)-k_dq_s$$

The dual-porosity model describes a macropore–micropore–solid interaction: macropores allow fast diffusion, while micropores slow dye migration due to adsorption and limited interphase exchange.

Simulation and Result

The coupled 3D diffusion–adsorption equations are solved numerically using FEM. The membrane is meshed into cubes subdivided into tetrahedra, with $C(\mathbf{x},t)$ approximated via linear basis functions in each element.

Figure 6 | FEM Algorithm Diagram

Time integration uses an implicit backward Euler scheme with Picard iterations to handle the nonlinear diffusion coefficient, updating solute concentration and adsorption at each step as: $$\mathbf{C}^{\,n+1} \approx (\mathbf{M} + \Delta t\, \mathbf{K}(\mathbf{C}^{\,n+1}_{(k)}))^{-1}\big(\mathbf{M}\mathbf{C}^n - \Delta t\,\mathbf{F}(\mathbf{q}^{\,n+1}_{(k)})\big)$$ $$\mathbf{q}^{\,n+1} = \frac{\mathbf{q}^n + \Delta t\, k_a\, \mathbf{C}^{\,n+1} q_{\max}}{1 + \Delta t (k_a \mathbf{C}^{\,n+1} + k_d)}$$ For intuitive simulation, stochastic particle tracking is used: dye release points are modeled as particles undergoing Brownian motion: $$\mathbf{x}_p^{\,n+1} = \mathbf{x}_p^{\,n} + \sqrt{2D_0 \Delta t}\,\boldsymbol{\xi}, \quad \boldsymbol{\xi}\sim \mathcal{N}(0, \mathbf{I}_3)$$ Reflective boundaries keep particles within the membrane, and their positions are tracked over time for dynamic visualization of dye diffusion.

Initially, the membrane contains no dye ($C(\mathbf{x},0) = 0$, $q(\mathbf{x},0) = 0$). FEM and iterative solution yield concentration profiles, adsorption dynamics, and particle diffusion trajectories.

Particle trajectories are simulated using Brownian motion, with reflective boundaries set to ensure particles remain within the membrane at all times. The rainbow colors do not affect the physical calculations; they serve solely as a “time-lapse visualization” to aid in understanding how particles diffuse through space over time.

Figure 7 | Single Particle Trajectory

As shown in the figure below, Brownian motion is used to introduce stochastic perturbations to the particles. In the large-pore region (the upper part of the plane), particles exhibit lower adsorption probabilities and therefore remain in continuous motion, represented by blue trajectories. In contrast, in the small-pore region, the adsorption probability is significantly higher, leading to frequent adsorption events that appear as red areas in the simulation.

Over time, as the simulation progresses, particles from the large-pore region gradually diffuse into the small-pore region and become adsorbed. Meanwhile, a certain probability of desorption also exists, allowing some particles to temporarily leave the adsorption sites. Ultimately, the system evolves such that most particles become adsorbed within the small-pore region, while a smaller fraction remains adsorbed in the large-pore region or continues to move within the small pores.

Figure 8 | BC's Dual-Porosity Simulation

By counting the number of particles, the 3D dye distribution is visualized through particle trajectory projection, showing rapid macropore transport and local delays in adsorption and slow interface exchange in micropores.

Figure 9 | 3D Dye Distribution

Particle-based results are linked to FEM simulations to evaluate dye penetration and adsorption kinetics. The model estimates saturation in ~37,825 s (~12 h), matching experimental data, supporting the FEM–particle approach for capturing macroscopic diffusion and microscale adsorption in BC membranes.

Figure 10 | Dyeing Simulation Result

• [1]Tan J, Zuniga C, Zengler K. Unraveling interactions in microbial communities - from co-cultures to microbiomes. Journal of Microbiology, 2015, 53: 295-305.
• [2]Ma L, Gao J P, Li D, Lian W Y. Dynamics of a delayed Lotka–Volterra competition model with directed dispersal. Real World Applications, 2023, 71: 103830.
• [3]Zhang Y, et al. Dyeing behavior of bacterial cellulose with natural dye coffee. PLoS ONE, 2022, 17(3): e0265743.
• [4]Morimoto K, et al. Computer simulation of dyeing process in textile design. Kyushu University, Design Faculty Reports, 2007.
• [5]Crank J. The Mathematics of Diffusion (2nd ed.). Oxford University Press, 1975.
• [6]Morimoto Y, Tanaka M, Tsuruno R, et al. Visualization of dyeing based on diffusion and adsorption theories. In: 15th Pacific Conference on Computer Graphics and Applications (PG'07), 2007, IEEE: 57-64.
• [7]Kim H, Kim H R. Production of coffee-dyed bacterial cellulose as a bio-leather and using it as a dye adsorbent. PLoS ONE, 2022, 17(3): e0265743.
• [8]Dhatt G, Lefrançois E, Touzot G. Finite element method. John Wiley & Sons, 2012.