Team:USTC

Overview

    To provide strategies for optimizing experimental results and to further verify the feasibility of the experimental design, we use the model as a prediction. In this project, the model we designed is divided into optimal enzyme ratio and diffusion simulation.

Optimal enzyme ratio

Introduction

    According to the design, we introduced a reaction pathway into the peroxisome of the Yarrowia lipolytica to produce borneol. It is a multi-step reaction from acetyl coenzyme A to GPP, the precursor of borneol. The relative amount of enzymes, i.e., enzyme ratio, plays a critical role in driving the overall reaction. In order to maximize the yield of borneol to provide adequate raw material for the drug, here we build a kinetic model combining five reactions from acetyl coenzyme A to mevalonate diphosphate to predict the optimal enzyme ratio for the reactions.

Description

    Let's start with the reaction pathway, cf. Fig. 1. In the peroxisome:

    The Yarrowia lipolytica produces acetyl coenzyme A through the β-oxidation of fatty acids in the peroxisome. Acetyl coenzyme A is then used to produce the precursor of borneol in an eight-step reaction that we have introduced. It is known that there are not above reactions and corresponding side reactions in the peroxisome before engineering. Theoretically, increasing the yield of mevalonate diphosphate would facilitate the production of GPP. The yield of mevalonate diphosphate is strongly linked to the enzyme ratio in the upstream reactions. So, we will explore the optimal enzyme ratio for the production of mevalonate diphosphate in the following section.

    First of all, we use Michaelis-menten equation to describe the kinetics of an enzymatic reaction:

    The following symbols are used in our models:



    Reaction of the synthesis of S6 and its upstream reactions:

    Reaction1:(E1)2S1 ⟺ S2+ CoA

    forward reaction:



    reverse reaction:



    Reaction2:(E2)S2 + CoA ⟹ S3



    Reaction3:(E3)S3 + 2NADPH + 2H+ ⟹ S4 + 2NADP+ + CoA



    Reaction4:(E4) S4 + ATP ⟹ S5 + ADP



    Reaction5:(E5) S5 + ATP ⟹ S6 + ADP



    Reaction of S6 consumption:

    Reaction6:(E6)S6+ ATP ⟹ S7 + Pi + CO2 + ADP



    The above parameters can be found in subsection Parameters.

    We ignore the reverse reactions of reactions 2 to 6 because their rate constants (k) are much smaller than forward reactions’. Using these rate equations, we can establish the relationship between the concentration of acetyl coenzyme A and the concentration of mevalonate diphosphate at certain enzyme ratio.

    If we want to explore the optimal enzyme ratio for reactions, we need to assume a constant value for the sum of enzymes and then change the enzyme ratio to compare the yield of GPP at different ratios. Based on the Genetic Algorithm, we search for the optimal resolution by simulating the natural evolutionary process.

    First, we set each enzyme to an initial value of 20 and fix the total number of enzymes to 120. We consider the relative content of the six enzymes as six genes and then code these genes into yeast ‘s “chromosome”. Each individual represents a solution for given problem.

    Second, we assume that a yeast produces one yeast as a generation and that Simple Mutation occurs in every generation, which means that a variation operation is performed on the values of two randomly specified motifs in the coding string, i.e. the values of two enzymes are randomly changed while the total number of enzymes remains unchanged.

    Third, for each offspring, we compare fitness with the previous generation and select. We decode the individual coding strings to obtain the phenotype, i.e. the enzyme ratio, and then we can calculate the corresponding S6 yield and compare the fitness of this generation with that of the previous generation. If it is higher than the fitness of the previous generation, we eliminate the previous generation and produce offspring based on this generation, otherwise we eliminate this generation and continue to produce offspring based on the previous generation.

    Thus each successive generation is more suited for their environment, which means that the optimal enzyme ratio can be calculated eventually.

    However, during the process, we found that if we used S6 yield directly as a function of fitness to compare yeast fitness, the reaction time steps would have to be within 1*10^ (-4) s in order to accurately obtain the concentration of the substrate in each reaction. Within a finite number of cycles, the program could only obtain S6 yield in a very short time, which couldn’t meet the yeast fermentation time designed.

    So, we made some improvements and used mathematical proofs to convert the problem of comparing S6 concentrations at equilibrium into a problem of comparing the slope of the S6 concentration as a function of time, so that the program could predict the S6 concentration at equilibrium in a shorter time.

    The following is the process of our proof.



    S6 concentration derivative with respect to time:




    When the concentration of S6 → +∞:



    Since the S6 concentration actually cannot tend to infinity, so:



    Solving the differential equation gives:



    When t → +∞, the reaction reaches equilibrium and the concentration of S6 tends to a constant value, from which we can deduce the concentration of S6 at equilibrium:



    Next, we prove that the greater the slope of the function of S6 concentration as a function of time, the greater the S6 equilibrium concentration.

    Clearly, the S6 equilibrium concentration is monotonically increasing with respect to C1.

    We consider S6 ,i.e. the slope at a concentration of S6 equal to half of the equilibrium concentration, when the slope has stabilised.




    Derivation with respect to C1 gives:



    Since the magnitude of C2 is around 10^(-6), we get:




    So, the slope of the function of S6 concentration as a function of time is monotonically increasing with respect to C1.

    Finally, we can conclude that the greater the slope of the function of S6 concentration as a function of time, the greater the equilibrium concentration of S6.

    Based on this, our program can predict and compare the magnitude of the S6 equilibrium concentration by the change in S6 concentration over a shorter period of time, and thus get the optimum enzyme ratio more efficiently.

Results


Fig.5: Optimal enzyme ratio

Parameters:

    Here we are going to introduce our choice of parameters.

    The following are the parameters that will be used in the rate equations.




    Considering that the carbon source of lipolytic yeast in the design is fatty acids and that the yeast peroxisome uses fatty acids for β-oxidation to produce relatively high concentration of acetyl coenzyme A, we considered the concentration of acetyl coenzyme A in the peroxisome as a constant value. The peroxisome is a relatively closed system and the major intermediates, end products and enzymes of our introduced reaction pathway were not present before the modification, so the initial S2 to S6 concentrations were 0. The enzyme concentrations for each reaction were on the order of μM.

Diffusion model

Introduction

    According to our design, we converted borneol to bornylamine and attached it to lipid nanoparticles and wanted to make it possible to cross the blood-brain barrier to exert drug action. To verify the feasibility of the design, here we model the diffusion of the lipid nanoparticles in the human body and predict the proportion of the lipid nanoparticles that would cross the blood-brain barrier.

Description

    Based on the multi-compartment model, we assume that the body consists of three parts: centralroom, peripheral room and blood-brain barrier. The central room includes blood and tissues that can be instantaneously distributed in balance with blood such as kidney, heart and liver. Peripheral room includes tissues with less blood supply such as fat, muscle, bone and cartilage.After administration, lipid nanoparticles are distributed to the central compartment immediately and then slowly diffuse to the peripheral room and bl-ood-brain barrier.


The following symbols are used in our models:

The following is our multi-compartment model, c. f. Fig1.

Fig.1: Multi-compartment model

    We assume that lipid nanoparticles enter the body in the amount of D0. Lipid nanoparticles cometo the central room, part of them are metabolized and excreted at a rate of v1, part of them enter the peripheral room at a rate of v12, and part of them come near the blood-brain barrier at a rate of v31. Lipid nanoparticles entering the peripheral room can likewise diffuse into central room at a rate of v21. A portion of lipid nanoparticles near the blood-brain barrier will diffuse into the central room at a rate of v31, and a portion will cross the blood-brain barrier at a rate of v3.


    Based on our model design, we can obtain the concentration of lipid nanoparticles at each site versus time and ultimately predict the proportion of lipid nanoparticles crossing the blood-brain barrier.

Results

    We assume that 1.5 g of lipid nanoparticles are administered intravenously into the body. We canget the relation between the concentration of lipid nanoparticles in each part versus time, cf. Fig.2.


Fig.2 Modeling image of the relation between concentration of lipid nanoparticles in each part and time

    We can predict that the concentration of liposomes crossing the blood-brain barrier is 3.8385·10^(-3)g/L , cf. Fig. 3.


Fig.3 Modeling image of the relation between concentration of lipid nanoparticles crossing the blood-brain barrier and time

Parameters

    Here we are going to introduce our choice of parameters. The following are the parameters that will be used in the diffusion model.


    Using data proven , we set K3=7.5·10^(-4)/min

References

[1] Hedl M, Sutherlin A, Wilding E I, et al. Enterococcus faecalis acetoacetyl-coenzyme A thiolase/3-hydroxy-3-methylglutaryl-coenzyme A reductase, a dual-function protein of isopentenyl diphosphate biosynthesis[J]. Journal of bacteriology, 2002, 184(8): 2116-2122.

[2] Middleton B. The kinetic mechanism of 3-hydroxy-3-methylglutaryl-coenzyme A synthase from baker's yeast[J]. Biochemical Journal, 1972, 126(1): 35-47.

[3] Primak Y A, Du M, Miller M C, et al. Characterization of a feedback-resistant mevalonate kinase from the archaeon Methanosarcina mazei[J]. Applied and environmental microbiology, 2011, 77(21): 7772-7778.

[4] Garcia D E, Keasling J D. Kinetics of phosphomevalonate kinase from Saccharomyces cerevisiae[J]. PLoS One, 2014, 9(1): e87112.

[5] Krepkiy D, Miziorko H M. Identification of active site residues in mevalonate diphosphate decarboxylase: implications for a family of phosphotransferases[J]. Protein Science, 2004, 13(7): 1875-1881.

[6] Verscheijden, L. F. M., Koenderink, J. B., de Wildt, S. N. & Russel, F. G. M. Development of a physiologically-based pharmacokinetic pediatric brain model for prediction of cerebrospinal fluid drug concentrations and the influence of meningitis. PLoS Comput Biol 15, e1007117, doi:10.1371/journal.pcbi.1007117 (2019).

[7] Lv, L. et al. Enhanced Anti-Glioma Efficacy by Borneol Combined With CGKRK-Modified Paclitaxel Self-Assembled Redox-Sensitive Nanoparticles. Front Pharmacol 11, 558, doi:10.3389/fphar.2020.00558 (2020).

ZJU-China ZJU-China