Highlights

1. Predict how the melanin production varies with time at different tyrosine (Tyr) concentrations.
2. Evaluate whether knocking out tyrA is suitable.

Background

    In melanin production process, Tyr is catalyzed by tyrosinase to synthesize dopaquinone. After a series of oxidation reactions, melanin can be produced. In this model, the series of dopaquinone oxidation reactions are simplified into a first-order chemical reaction. Kinetics are used to establish the reaction mechanism to describe the production and consumption of chemicals, allowing us to predict how the melanin production changes over time.

    To prevent bacteria from leaking out, tyrA gene, which is essential for Tyr synthesis in E. coli, is knocked out, so the bacteria can only take in Tyr from the medium and are unable to live outside (see Safety page). To confirm that deletion of tyrA will not inhibit the bacteria growth and hinder the melanin production, three different Tyr concentrations were chosen to simulate the consumption curves in our model.

    In general, we eliminated the insignificant parameters, focused on only the main reactant, Tyr, and successfully predicted the melanin yield by modeling.

Melanin Pathway Model

I. Reaction mechanism

    Melainin pathway equation:

\(\ce{T+E <=>[k_1][k_2] ET ->[k_3] QH+E}\)    \(\ce{QH ->[k_4] M}\)

Table 1. Descriptions of melanin pathway parameters

Parameter	Description	Unit
\[[E]\]	Tyrosinase concentration	\[g \cdot L^{-1}\]
\[[T]\]	Tyr concentration	\[g \cdot L^{-1}\]
\[[ET]\]	The concentration of the combination form of Tyr and tyrosinase	\[g \cdot L^{-1}\]
\[[QH]\]	Dopaquinone concentration	\[g \cdot L^{-1}\]
\[[M]\]	Melanin concentration	\[g \cdot L^{-1}\]
\[k_1\]	Rate constant for melanin pathway equations	\[L \cdot g^{-1} \cdot hr^{-1}\]
\[k_2, k_3, k_4\]	Rate constants for melanin pathway equations	\[hr^{-1}\]

Table 2. Descriptions of melanin pathway equations

Equation	Description	Unit
\[{d[T] \over dt} = {-k_1[E][T]+k_2[ET]}\]	The consumption rate of Tyr	\[g \cdot L^{-1} \cdot hr^{-1}\]
\[{d[QH] \over dt} = {k_3[ET]}-k_4[QH]\]	The formation rate of dopaquinone	\[g \cdot L^{-1} \cdot hr^{-1}\]
\[{d[E] \over dt} = {-k_1[E][T]+k_2[ET]+k_3[ET]}\]	The rate of changing tyrosinase concentration	\[g \cdot L^{-1} \cdot hr^{-1}\]
\[{d[ET] \over dt} = {-k_3[ET]}+k_1[E][T]-k_2[ET]\]	The formation rate of the enzyme-substrate combination form	\[g \cdot L^{-1} \cdot hr^{-1}\]
\[{d[M] \over dt} = {k_4[QH]}\]	The production rate of melanin	\[g \cdot L^{-1} \cdot hr^{-1}\]

II. Assumptions

1. The initial Tyr concentration in the medium is equal to the initial intracellular Tyr concentration.
2. The series of dopaquinone oxidation processes is negligible and was viewed as a first-order chemical reaction.

III. Results

    In Fig. 1, after 5 days of culture, there were still 10% of Tyr left. This indicated that 0.05 g/L of Tyr is sufficient for melanin production during cell culture. In conclusion, melanin synthesis will not lead to Tyr depletion with the supply of at least 0.05 g/L of Tyr. Moreover, extra Tyr can support the growth of bacteria with tyrA deleted, which cannot synthesize Tyr itself.

Fig. 1. Melanin production curve (adding 0.05 g/L Tyr)

Fig. 2. Melanin production curve (adding 0.1 g/L Tyr)

Fig. 3. Melanin production curve (adding 0.2 g/L Tyr)

     In Table 3, unlike chemical reactions, the rate constants are different when Tyr concentration varies. In addition, a confirmation test was conducted to ensure the external validity of the model (see Appendix).

Table 3. Estimated parameter values in melanin pathway

Rate constant	0.05 g/L Tyr	0.1 g/L Tyr	0.2 g/L Tyr
\[k_1\]	2.506	7.576	11.443
\[k_2\]	0.576	5.647	6.049
\[k_3\]	0.144	0.134	0.209
\[k_4\]	0.029	0.041	0.030

Conclusions

1. Melanin production can be predicted.
2. Melanin synthesis will not lead to Tyr depletion with the supply of at least 0.05 g/L of Tyr. Moreover, extra Tyr can support the growth of bacteria with tyrA deleted.

    To improve our model in the future, we hope to measure melanin yield more precisely by melanin extraction to better simulate the production curve.

Highlight

    Taguchi methods are used to optimize the experimental parameters and successfully increased melanin yield by 291%.

Background

    In addition to Tyr concentrations, many other factors may also significantly affect melanin yield. However, if too many factors were included in the kinetic equation, it would take much time and effort to find multiple kinetic parameters and build models. Taguchi methods, developed by Genichi Taguchi, are a popular quality control approach in manufacturing used to produce high-quality products and minimize quality variations. Taguchi methods use orthogonal arrays to maximize test coverage by pairing and combining the inputs and test the system with less number of cases for time-saving. This is a powerful approach to optimize multiple factors at a time and minimize the number of experiments.

Why Choose Taguchi Methods?

Table 4. Comparison of the common quality control methods, including traditional experimental design, factorial design, RSM, and Taguchi methods. Factorial design, RSM, and Taguchi methods are all popularly used in design of experiment (DOE).

Methods	Our consideration
Traditional experimental design	Lack efficiency due to too many experiments Limited numbers of control factors and levels Continuous variables are required in RSM. Compared to Taguchi methods, factorial design and RSM are much more complicated to understand.
Factorial design
Responsive surface methodology (RSM)
Taguchi methods	Variables can be continuous or discrete. Allows more control factors and levels within minimum number of experiments Simple and efficient

    Taguchi methods allow simultaneous examination of multiple factors and extract quantitative information by conducting only a few experiments.

Taguchi Methods

I. Experiment design

Fig. 4. Taguchi methods step 1-4.

Step 1. Selection of quality characteristics, control factors, and levels

1. Quality characteristic: OD₄₀₀ value (OD₄₀₀ is the optical density at wavelength 400 nm. One OD₄₀₀ unit is equivalent to 0.066 g/L of eumelanin[2]).

2. Table 5. Control factors and levels

   	A: Tyr addition time	B: Cu2+ conc. (mM)	C: LB addition (%)
   Lv 1	OD600 = 0.6	0.2	0
   Lv 2	day 1	0.5	1
   Lv 3	day 2	0.8	10

3. Why did we choose the above control factors and levels?
   • Tyr addition time
         E. coli does not take in much Tyr. Bacteria in different growth phases probably consumed Tyr to different extents. Hence, we considered the Tyr addition time as one of the control factors.
   • Cu²⁺ conc. (mM)
         Tyrosinase is a copper-containing enzyme. Low Cu²⁺ concentration promotes tyrosinase activity, while high copper concentration is lethal to cells. Experiment results showed that 0.5 mM could yield the most melanin (see Results page). Therefore, we selected 0.2 mM, 0.5 mM, and 0.8 mM.
   • LB addition in M9 medium (%)
         M9 has a higher C/N ratio and is commonly used for melanin synthesis in several previous studies, while LB medium has a lower C/N ratio and is more nutritious for E. coli (see Results page). We decided to add 0%, 1%, and 10% of LB into M9 medium.

Step 2. Orthogonal design of experiments

    We selected 3 levels each. Using Taguchi L9 orthogonal array, only 9 experiments were required to be conducted. The orthogonal array helped us strategically design test cases to balance out the effect of each factor.

Table 6 & 7. L9 orthogonal array (left); L9 orthogonal array for experimental design of melanin synthesis (right)

A B C L1 1 1 1 L2 1 2 2 L3 1 3 3 L4 2 1 2 L5 2 2 3 L6 2 3 1 L7 3 1 3 L8 3 2 1 L9 3 3 2

A: Tyr addition time B: Cu2+ conc. (mM) C: LB addition (%) L1 OD600 = 0.6 0.2 0 L2 OD600 = 0.6 0.5 1 L3 OD600 = 0.6 0.8 10 L4 day 1 0.2 1 L5 day 1 0.5 10 L6 day 1 0.8 0 L7 day 2 0.2 10 L8 day 2 0.5 0 L9 day 2 0.8 1

Step 3. Taguchi experiment & Data analysis

    Experiments were conducted using the Taguchi orthogonal design. The result data was analyzed to determine the optimal factor level.

Step 4. Confirmation experiment

    A confirmation test is required for the following two reasons: First, to test the reproducibility of Taguchi experiments. Second, to prove that our optimized parameter can indeed increase melanin yield. Therefore, the confirmation test is designed with two experiments: any experiment picked from L1 to L9 (L6 was picked in our test.) and one with the optimized parameter.

II. Results

Table 8. OD₄₀₀ results for L1 to L9 experiments

	A: Tyr addition time	B: Cu2+ conc. (mM)	C: LB addition (%)	OD400
L1	OD600 = 0.6	0.2	0	4.102
L2	OD600 = 0.6	0.5	1	5.155
L3	OD600 = 0.6	0.8	10	4.876
L4	day 1	0.2	1	5.594
L5	day 1	0.5	10	5.168
L6	day 1	0.8	0	2.744
L7	day 2	0.2	10	4.273
L8	day 2	0.5	0	3.504
L9	day 2	0.8	1	1.672

Table 9. Response table for means

	A: Tyr addition time	B: Cu2+ conc. (mM)	C: LB addition (%)
1	4.711	4.656	3.450
2	4.502	4.609	4.141
3	3.150	3.098	4.772
Delta (= Max - Min)	1.561 (Lv1 - Lv3)	1.559 (Lv1 - Lv3)	1.322 (Lv3 - Lv1)
Significance Rank	1	2	3

• Mean value: Take Cu²⁺ 0.2 mM for example. L1, L4, L7 experiments all contain this parameter. The OD₄₀₀ of L1 (4.102), L4 (5.594), L7 (4.273) were averaged to get the mean value 4.656 for Cu²⁺ 0.2 mM.
• Delta value: Difference between the highest and the lowest characteristic average of a particular factor. A larger delta value indicates that the level difference of the factor has a greater effect on the quality characteristic, which means the control factor is more significant.

    In Fig. 5, adding Tyr when OD₆₀₀ = 0.6 (A-L1), 0.2 mM Cu²⁺ (B-L1), 10% LB (C-L3) is the optimal condition to increase melanin production. Compared to LB addition, Tyr addition time and Cu²⁺ conc. are more signifcant in the Taguchi experiment.

     For the final step, a confirmation experiment was conducted to verify the reproducibility of Taguchi experiments and the melanin yield under the optimized conditions.

Fig. 5. Response graph for means

    In Table 10, the OD₄₀₀ value from the confirmation test is close to that in the Taguchi experiment. The result verified the reproducibility of our Taguchi methods. In addition, the experiment with optimized parameters had a higher OD₄₀₀ value than all the Taguchi experiments.

Table 10. Confirmation experiment result

	A: Tyr addition time	B: Cu2+ conc. (mM)	C: LB addition (%)	OD400	Taguchi L6 OD400
L6	day 1	0.8	0	2.824	2.744
Optimized	OD600 = 0.6	0.2	10	5.940

Conclusions

    In Table 11, we compared the melanin yield under optimized condition (0.239 g/L) and original condition (0.061 g/L). It proves that by using Taguchi approach, we successfully increased melanin yield from previous wet lab experiments by 291%. Moreover, our wet team conducted experiments to compare the melanin production from dual plasmids under the optimized and original condition. In Table 12, under the optimized condition, we also increased our melanin yield significantly.

Table 11. Comparison of melanin yield produced in different culture conditions

	Medium	Temp. (℃)	Cu2+ conc. (mM)	Host	Tyr conc. (g/L)	Gene	Time (hr)	Melanin (g/L)
Original	M9	37	0.5	DH5α	0.4	pSUI-Ptrc-melA-B0015	120	0.061
Optimized	M9 + 10% LB	37	0.2	DH5α	0.4	pSUI-Ptrc-melA-B0015	121	0.239
							94	0.189
Reference[3]	M9	30	0.315	W3110MGT	0.4	PtrcMut-melA	90	0.050

Table 12. Comparison of melanin yield produced from dual-plasmid bacteria (melA and gadB)

	Medium	Temp. (℃)	Cu2+ conc. (mM)	Host	Tyr conc. (g/L)	Gene	Time (hr)	OD400
Original	M9	37	0.5	DH5α	0.4	Dual plasmids (melA and gadB)	82	6.190
Optimized	M9 + 10% LB	37	0.2	DH5α	0.4			8.214

Highlights

1. Establish mathematical relationships about how the bacteria accumulate with time as a consequence of the nutrient uptake.
2. The rate constants were obtained for both the bacterial accumulation and the nutrient uptake, allowing us to predict the time required to replace the medium.

Background and Motivation

    The observed bacterial accumulation was reasoned by us as a consequence of the nutrient uptake. Both processes are determined simultaneously by the instantaneous bacterial and nutrient concentrations. A kinetic model is necessary to quantify these processes, allowing us to find the rate constants for these processes from the experimental data. It also tells when to re-supply the nutrients that are consumed by the bacteria.

Kinetic Model for Nutrient Uptake and Bacterial Accumulation

    NCKU_Tainan proposed a simple kinetic model to express the bacterial accumulation rate and the nutrient uptake rate in terms of the bacterial concentration and the nutrient concentration. The unknown rate constants here were obtained by fitting the bacterial accumulation curve measured from the experiment to the one predicted by the model. This model was used to compute the nutrient uptake rate, allowing us to find the suitable medium exchange time for sustaining the bacterial accumulation. In addition, several studies have shown that 3D cell structure has a larger surface area resulting in higher nutrient uptake rate, which is expected to occur in the hanging drop microfluidic chip[5].

Table 13. Descriptions of parameters for nutrient uptake and bacterial accumulation rate calculation

Parmeter	Description	Unit
\[u\]	Bacteria concentration	\[g \cdot L^{-1}\]
\[v\]	Nutrient concentration	\[g \cdot L^{-1}\]
\[k\]	Rate constant of bacterial accumulation	\[L \cdot g^{-1} \cdot s^{-1}\]
\[R \times k\]	Rate constant of nutrient uptake The rate of nutrient uptake is k times larger than the bacterial accumulation rate.	\[L \cdot g^{-1} \cdot s^{-1}\]

Table 14. Descriptions of equations for nutrient uptake rate and bacteria accumulation rate

Equation	Description	Unit
\[{du \over dt} = {kuv}\]	Bacterial accumulation rate Bacterial accumulation is assumed to be positively correlated to nutrient uptake.	\[g \cdot L^{-1} \cdot s^{-1}\]
\[{dv \over dt} = {-Rkuv}\]	Nutrient uptake rate Since a single bacterium consumes more nutrients than how much it needs to grow properly, nutrient uptake rate should be faster than the bacteria accumulation rate.	\[g \cdot L^{-1} \cdot s^{-1}\]

I. Assumptions

1. The bacterial accumulation rate is proportional to the bacterial concentration and the nutrient concentration.
2. The same proportionality as the above is also applied to the nutrient uptake rate but with a larger rate constant.

II. The Model

    Using the assumptions above, a coupled set of differential equations were constructed for the bacterial concentration (u) and the nutrient concentration (v).

where

\[{Q}={u_0 + {v_0 \over R}}\text{ ......(Eq. 3)}\]

with u₀ and v₀ being the respective bulk concentrations for the bacteria and the nutrients.

    Eq. 1 indicates that the bacterial concentration u approaches toward Q as t~∞, allowing us to find the value of Q from the long-term bacterial concentration observed in our experiment. Also because both u₀ and v₀ are known, the value of R can also be determined and the remaining rate constant k can be found by fitting the observed bacterial accumulation data to Eq. 1. After Q, R, and k values are obtained, nutrient uptake with time can be modeled by using Eq. 2.

III. Results

    Table 15 lists the fitted values of the rate constants k and Rk from the culture experiment in the microfluidic chip. The results are also compared to those in the no drop formation case. These constants in the chip case are higher than those in the no drop formation case, indicating that the microfluidic chip did expedite the bacterial accumulation and nutrient uptake (see Appendix).

Table 15. Comparison between the values of k and Rk obtained from our microfluidic chip and those in the no drop formation case

	\[k\text{ }(L \cdot g^{-1} \cdot hr^{-1})\]	\[R \times k\text{ }(L \cdot g^{-1} \cdot hr^{-1})\]
Chip without hydrophilic membrane	0.030	6.5
Chip with hydrophilic membrane (No drop formation)	0.017	6.0

    Fig. 6 displays the growing bacterial concentration data collected from our microfluidic chip, together with the fitted model curve for capturing the data trend. Fig. 7 plots the nutrient uptake curve predicted by the model. It shows that after 8 hours of the bacterial culture the nutrient concentration level is reduced to 3.2% of the supplied nutrient concentration (0.48 g/L), suggesting the need to re-supply the medium (see Hardware page).

Fig. 6. Bacterial concentration data collected from the microfluidic chip and the fitted model curve

Fig. 7. Temporal evolution of the nutrient concentration predicted by the model

Conclusions

1. The proposed kinetic model can successfully capture the growing data trend of the measured bacterial concentration.
2. By fitting the experimental data with the model’s predictions, the rate constants for both bacterial accumulation and nutrient uptake can be found. This also allows us to demonstrate the facilitation of these two processes by the patterned microfluidic chip.

Highlights

1. The modified Keller-Segel model is capable of providing a quantitative account for the spatiotemporal behavior of the bacterial growth in an incubation microwell in the microfluidic chip.
2. The model was able to predict the extent of the bacterial growth after a few hours of incubation period, successfully capturing the experimental observation after 8 hr of bacterial culture in our microfluidic chip.

Background and Motivation

    The kinetic model has demonstrated that the bacterial accumulation can be promoted with the aid of the microwells in the microfluidic chip (Fig. 8). If we can know how the bacteria grow in both space and time locally within a microwell, this will provide an added advantage of optimizing the spatiotemporal growth of the bacteria via changing the size of the microwell.

Fig. 8. Apparent bacterial aggregation (white spots) in the microwells after 8 hr of bacterial culture

    To describe the spatiotemporal growth of the bacteria, it is necessary to use a transport model for capturing the spatial distributions of the bacteria and the nutrients and how these distributions vary with time. Since any concentration gradient naturally causes diffusion from high to low concentration regions, this effect will oppose the bacteria buildup and the nutrient depletion. Therefore, aside from the inherent kinetics-driven bacterial accumulation and nutrient uptake, the model has to further include diffusion for both the bacteria and the nutrients so that a steady state can be established due to gradual slowdown of these accumulation and consumption processes.

    In fact, such a diffusion-mediated accumulation/consumption effect is quite like the one in chemotaxis in which microorganisms move toward food-rich (or run away from toxin-rich) places. Knowing that the chemotaxis system has been successfully described by the Keller-Segel model[6], we modified this model to provide a quantitative account for the observed bacterial growth in an incubation microwell in our microfluidic chip.

Bacteria Aggregation Simulating Model

I. Modified Keller-Segel Equation

    By simplifying the equation, we could better bridge the gap between the simulation results and experimental observations.

    At the beginning of bacteria culture, bacteria are unevenly distributed. Around the small bacteria clusters, nutrients are consumed faster, and the bacteria can thus grow faster. As time goes by, bacteria aggregation can be observed in a macrocosm vision.

Table 16. Parameter descriptions and values

Parameter	Description	Value
\[D_1\]	Diffussion coefficient of bacteria	1.5 × 10-12 m2/s[7]
\[D_2\]	Diffussion coefficient of nutrients	1.5 × 10-9 m2/s[7]
\[v_0\]	Initial nutrient concentration	15 g/L
\[u_0\]	Initial bacteria concentration	0.269 g/L
\[L\]	Radius of bacterial aggregates	200 μm
\[x_0\]	Well radius	0.5 mm
\[k\]	Rate constant of bacterial growth	0.00005064 L ⋅ g-1 ⋅ s-1
\[t^*\]	The characteristic time	7.4 hr
\[kuv\]	Bacterial accumulation rate
\[Rkuv\]	Nutrient uptake rate

II. Assumptions

1. The bacterial growth and aggregation behavior in one microwell can well represent the average bacterial behavior in all microwells in the microfluidic chip.
2. Use the bacterial concentration to describe the extent of the local bacterial buildup and crowding.
3. The local bacterial buildup caused by the nutrient uptake is counteracted by the natural tendency to lower the bacterial concentration due to the bacteria’s self-diffusion.
4. The nutrient depletion around a bacterial aggregate is also compensated by the nutrients' diffusion from the bulk to the aggregate.

III. Simulation method

    The simulation was conducted by using MATLAB for numerical integration of the modified Keller-Segel equations (see Appendix). The equations, variables and parameters were expressed in the dimensionless form (see below and Table 17). This provided a more convenient way to measure variations of physical quantities with respect to their characteristic values. This also helped us to identify relative importance between relevant effects.

    After obtained from the simulation, the results were converted back to dimensional ones by multiplying them by their characteristic values. e.g. dimensionless time = 1, characteristic time = 7.4 hr, thus, the actual time is 1 × 7.4 hr.

    Dimensionless expression:

Table 17. Dimentionless form of variables

Variable	Description	Unit	Dimensionless	Characteristic quantity
\[u\]	Bacteria concentration	\[g \cdot L^{-1}\]	\[U\]	\[{U}={u \over u_0}\]
\[v\]	Nutrient concentration	\[g \cdot L^{-1}\]	\[V\]	\[{V}={v \over v_0}\]
\[t\]	Time (aggregation time)	\[hr\]	\[t^*\]	\[{t^*}={L^2 \over D_1}\]
\[x\]	Distance	\[mm\]	\[X\]	\[{X}={x \over x_0}\]
\[k\]	Rate constant of bacterial growth	\[L \cdot g^{-1} \cdot hr^{-1}\]	\[K\]	\[{K}={k{L^2 \over D_1}v_0}\]
\[t\]	Time	\[hr\]	\[T\]	\[{T}={t \over t^*}\]
\[D_1\text{, }D_2\]	Diffusion coefficients of nutrients and bacteria	\[m^2 \cdot s^{-1}\]	\[D\]	\[{D}={D_2 \over D_1}\]

IV. Results

    MATLAB was used to convert the solution into dimensional modified KS equation and plot a 3D picture. In Fig. 9, bacteria concentration increased over time and decreased over the distance from the origin.

Fig. 9. Distribution of bacteria in a well

    As shown in Fig. 10, bacteria concentration is the highest at the center of the well where bacteria aggregation takes place (distance = 0). On the other hand, gradual decrease of the bacteria concentration occurred as it gets closer to the edge.

    In Fig. 11, the blue area was defined as the bacteria aggregation area, and the red line represents the radius of the well. The diameter of the aggregation area is approximately 40% of the well diameter, which is 200 μm. The highest bacteria concentration value is 0.11, 0.11 - (0.11×10%) = 0.1 in Fig. 9. The area with bacteria concentration above 0.11 occupied of 36% of well diameter. In conclusion, the simulation result corresponds to our experimental observations.

Fig. 10. Bacteria distribution in a well after 2, 4, 6, and 8 hr (left)
Fig. 11. Bacteria aggregation in the microfluidic chip after culturing 8 hours. Paint 3D was used, and the color tolerance was set to 10% to select pixels with darker colors. (right)

    Through dimensional analysis, we could calculate the characteristic time t* = L² / D₁ = 200² / (1.5 × 3600) = 7.4 hr. ( L is the radius of bacteria cluster, and D₁ is the diffusion coefficient.)

    In Fig. 12, bacteria concentration has come to steady state when T = 1, which corresponds to 7.4 hr. In our hardware experiment, bacterial aggregation could be easily observed in the well after around 6-8 hours (see Hardware page). This indicates that our model matches the experimental results.

Fig. 12. Bacteria concentration at the center of a well

Conclusions

    Transport Model for Bacterial Growth provides an insight into the formation of bacteria clusters inside the well and described bacteria aggregation in a quantitative way.

    For further studies, there is a possibility of capturing microbial motion, providing animations to demonstrate how bacteria move and aggregate. This surely is a promising way for model modification and advancement.

    Studies have shown that the aggregation mechanism can make bacteria much less susceptible to antibiotics in space; thus, astronauts have a higher chance to be infected with serious bacterial diseases in space than on Earth[8, 9]. By using this model built by NCKU_Tainan, bacteria aggregation in space can be predicted to determine when to sterilize the equipment in space stations before bacteria clusters grow larger and have higher antibiotic tolerance. Therefore, we can not only provide a solution to astronauts’ health issues but also let people safely travel to space in the future.