Modeling supports the use of tyrosinase over laccase for melanin production

First, to save time and allow for faster screening, a double substrate model was used to analyze two potential candidates for melanin production in Saccharomyces cerevisiae. Laccase and tyrosinase have both been identified as enzymes used by microbes to produce melanin from the amino acid tyrosine or its metabolites (Tran-Ly et al., 2020). One drawback of laccase is that it requires L-DOPA as a precursor for melanin. This limitation could be avoided by using tyrosinase. In this case, tyrosine, which in contrast to L-DOPA, is abundant in yeast cells, could be used as the starting molecule for melanin synthesis (Eisenman & Casadevall, 2012). However, we wanted to evaluate if it was beneficial to introduce both enzymes. In this case, tyrosinase would synthesize L-DOPA, a substrate for laccase that converts it to dopaquinone. To test the effectiveness of this strategy, we decided to use a double substrate ODE model to model a single oxidation reaction of L-DOPA to dopaquinone when catalyzed by either tyrosinase or laccase (Fig. 1).

In this model the two substrates, L-DOPA and oxygen, were the two limiting factors of the reaction. We kept the enzyme concentration constant in these models, because the aim was to compare the reaction rate when catalyzed by different enzymes.

The general equation for the double substrate model of enzyme-catalyzed L-DOPA oxidation kinetics was introduced:

As the basis for building the model we used the constant of maximal reaction rate per mass of enzyme (Vm, U mg^–1), Michaelis-Menten constants for L-DOPA and dissolved oxygen (Km ldopa and Km^O2, mmol dm^–3), constants that denote concentrations of L-DOPA (Cldopa, mmol dm^–3), dissolved oxygen (C_O2, mmol dm^–3) and enzyme (Yenzyme, mg cm^–3). The values for constants are presented in Supplementary Table 1.

For all our models, we used an approximation of tyrosinase enzyme concentration taken from literature because it is a realistic value for other organisms (Tahany M. Abdel-Rahman, 2019).

The model revealed two significantly different concentration-rate dependencies of L-DOPA oxidation (Fig. 1), which corresponded to the reaction being catalyzed by either laccase or tyrosinase. Tyrosinase has a substantial advantage in reaction velocity compared to laccase. We find that due to the marked difference in the reaction rates for the two enzymes, the contribution of laccase in a two-enzyme approach would be small. Hence, we decided to continue only with tyrosinase as the main enzyme of the melanin pathway.

Modeling shows the effect of oxygen concentration on the reaction rate

As mentioned in our ENGINEERING page, we tried different strategies to achieve the most efficient melanin production. After we carried out the EXPERIMENTS, we realized that extracellular melanin production is better than intracellular. We hypothesized that this could be because of a more reducing environment inside the cells, while melanin synthesis contains oxidation steps. So, we modified the previous model to test how oxygen concentration affects the L-DOPA oxidation reaction rate. The detailed description of the model mathematics is the same as for Figure 1. It is important to mention that the range of dissolved oxygen concentration value is set up according to the estimated parameters (from 0 to 1 mmol dm^–3) (Tišma et al., 2008).

The simulation shows a strong correlation of the reaction rate and oxygen concentration. A tendency of rapid increase in the reaction rate with the increase of dissolved oxygen concentration on each point of L-DOPA concentration can be observed (Fig. 2). This model supports our hypothesis based on the experimental data, and encourages us to consider the redox state of the environment as a key parameter in melanin synthesis.

Modeling to evaluate the potency of melanin production in yeast

Next, we moved on to estimate the levels of melanin we could reasonably expect from the strain. For this, we proceeded to model the entire eumelanin biosynthesis pathway. This introduced certain difficulties, as the earlier steps of the pathway are enzymatically catalyzed, while later steps are spontaneous and their rates are difficult to predict (Eisenman & Casadevall, 2012). For this reason, we modeled the reaction of the pathway up to dopaquinone synthesis using ODE’s. All reactions following dopaquinone were modeled using an Agent-Based (AB) model. Due to the stochastic nature of AB models, we expect them to have a better chance to predict the spontaneous steps of the pathway.

We used the Michaelis-Menten equation for an ODE model to estimate the pathway steps from tyrosine to L-DOPA and L-DOPA to dopaquinone (Fig. 3):

where for our model E is tyrosinase enzyme, S is tyrosine substrate, ES complex represents L-DOPA intermediate, P - dopaquinone.

We used SimBiology library in MATLAB to generate the model. Firstly, we introduced all compounds, such as tyrosine, L-DOPA, and dopaquinone, to the above equation. Initial values of each of them were defined as a concentration in mg/L. The two reactions with corresponding kinetic laws and parameters were defined as follows:

1) Tyrosine conversion into L-DOPA, driven enzymatically by tyrosinase. MassAction kinetic law and reaction rate constants k1 (k^on) and k1r (k^off) expressed in (mole×second)^–1 and s^–1 respectively. K1r value was not required, since it was automatically calculated by the toolkit.
2) Doapquinone synthesis from tyrosine, catalyzed by tyrosinase. MassAction kinetic law and reaction rate constant k2 (k^cat) was expressed in s^–1.

The values for these constants are shown in Supplementary Table 1

After defining all necessary features for the model, we obtained such rate rules (ODEs) for each species:

Using sbiosimulate function we built the graph which represents the concentration of each species in time (Fig. 3). The model indicated that L-DOPA reaches its maximal concentration of 16.2 mg/L in two minutes. It also estimates the amount of an important intermediate, dopaquinone, which is the last metabolite derived from the enzymatically-catalyzed reactions. This is an input for our subsequent AB model.

For the AB model the following pathway was introduced:

dopaquinone → cyclodopa → dopachrome → dhi → indolequinone → eumelanin

The AB model allows us to estimate product concentration using initial substrate concentration and reaction rates for each of the intermediate reactions.

The core principles and techniques for building AB models are uniform among biological modeling languages. We used the Kappa language, because it has relatively simple logic and syntax.

To build the AB model, we first introduced all agents (substrate with the states of dopaquinone, cyclodopa, dopachrome, dhi, indolequinone, eumelanin, and agents that represent each reaction), and constants, which will be described later in this section. For the initial concentration of dopaquinone, we used data from the ODE model (see the previous section). In the stochastic model, we are using countable quantities of agents. So we need to convert values of agents’ molar concentrations to the numbers of molecules in the cell. This requires the determination of the cell volume. For our simulation, we decided to model a fraction of the whole cell in order to save time and computing resources. A typical eukaryotic cell volume is about 2.5 x 10^-12 L.

We used Avogadro's Number and the model volume to convert molar concentrations to a number of molecules:

where n is the number of molecules of the agent a with molar concentration [a]. A is Avogadro’s Number, and V is a fractional cell volume that we used in our model.

There are several types of reaction modeling systems in the Kappa language (Fig. 4).

For every reaction in the pathway, we used several mechanisms, namely binding, unbinding, and modification. As these reactions are not catalyzed by enzymes and occur spontaneously with the presence of the substrate, binding/unbinding can be considered as close to zero. Even with this limitation, as it requires binding and unbinding constants, the model is perfectly suited to simulate the spontaneous steps of the pathway through stochastic modeling. The rate constants for “modification” mechanism were set according to the literature (Edge et al., 2006; Ito & Wakamatsu, 2008; Land et al., 2003; Odh et al., 1993).

Using this model, we built a graph that illustrates the dependence of the number of synthesized eumelanin molecules from time (Fig. 5A). For better visualization of the stochastic nature of the model, (Fig. 5B) shows a selected part of the reaction time course. The model predicts melanin production of around 14 000 molecules by one cell when production reaches its steady state. This is about 1.767 x 10^-5 mg/μL of melanin produced in a single cell. For comparison, a recent study showed production of around 0.420 g/L melanin in an optimized cell factory using black yeast (Elsayis et al., 2022). This modeling shows that while there is room for improvement in melanin production in S. cerevisiae compared to black yeast, we could expect reasonable production efficiency in our yeast strains.