Enzymes and metabolism

Life runs on chemistry

In Lesson 3 you followed a gene all the way to a protein. Most of the proteins a cell makes are enzymes: molecules that make one chemical reaction happen, over and over. This lesson is about what they do. By the end you will have written Python that predicts how much energy syn3A can make from its food, and which of its genes it cannot live without.

Metabolism is the name for all the chemical reactions in a cell. Some break food down to release energy. Others use that energy to build the cell’s parts: amino acids, nucleotides, lipids. Think of the cell as a factory. Raw materials come in through the doors, machines (the enzymes) turn them into parts, and the factory runs on one kind of power.

For syn3A, someone has written the whole factory down. In 2019 Marian Breuer, Zaida Luthey-Schulten and their colleagues published a metabolic model of syn3A called iMB155: a list of every reaction they could find evidence for. It has

• 338 reactions: 85 of them are “exchanges” that let a molecule enter or leave the medium around the cell, 86 carry molecules across the membrane, and 167 change one molecule into another inside the cell;
• 304 metabolites (small molecules that take part in reactions): 220 inside the cell and 84 outside;
• 155 genes, whose proteins run the reactions. That is about a third of syn3A’s protein genes.

Inside the cell, the reactions sort into a few areas: 50 for nucleotides, 43 for amino acids, 32 for central metabolism (breaking down sugar), 19 for lipids, 14 for vitamins and other helpers, 8 that assemble the big molecules, and 1 that adds up everything the cell needs to grow. This lesson uses the model’s real data, so you will see this credit wherever it appears: Metabolic model iMB155 from Breuer et al. 2019, eLife 8:e36842, CC BY 4.0.

What syn3A imports and what it makes

A bacterium like E. coli can build almost everything it needs from sugar and a few salts. syn3A can’t. Its genome was cut down to the bare minimum, and most of the genes for building blocks went with it. So it imports almost every part ready-made. In the model, 78 kinds of molecule can come in:

• glucose, its only source of energy, and glycerol, a building block for lipids;
• all 20 amino acids, loose, and most of them also as short chains of four;
• nucleosides and bases, the parts of DNA and RNA letters;
• fatty acids and cholesterol for the membrane;
• vitamins and pieces of helper molecules, such as riboflavin, nicotinate and coenzyme A;
• salts, ammonia, oxygen, carbon dioxide and water.

What does it send out? Mostly lactate and acetate, the leftovers of burning glucose, plus some pyruvate and carbon dioxide.

You may have learned that cells burn sugar all the way to carbon dioxide, using oxygen, in two big stages called the TCA cycle (or Krebs cycle) and the respiratory chain. That is how your own cells get most of their energy. syn3A has neither. It gets all of its energy from the first, short part of sugar breakdown, glycolysis, and throws the half-used leftovers away as lactate and acetate. As the 2022 whole-cell model put it, “all ATP in Syn3A is generated by the central metabolism.” It is wasteful, but it needs only about a dozen proteins. You will build that pathway, gene by gene, later in this lesson.

Enzymes make reactions fast

Many reactions in a cell are “downhill”: the products have less energy than the starting molecules, so the reaction can happen on its own. But downhill doesn’t mean fast. In plain water at room temperature, a peptide bond (the link between two amino acids in a protein) takes about 400 years to break by itself. Pulling a phosphate off a sugar takes about a million years. A cell that divides every two hours can’t wait that long.

A catalyst is a molecule that makes a reaction go faster and comes out unchanged at the end, ready to do it again. An enzyme is a catalyst made of protein. Enzymes often speed reactions up by ten orders of magnitude (10¹⁰ times) or more.

Each enzyme has a pocket, the active site, shaped to hold one particular molecule, its substrate. The enzyme turns the substrate into a product and lets it go. In the 1890s Emil Fischer pictured the pocket as a rigid lock with the substrate as its key. Today we know many enzymes move: they close around the substrate once it has docked, like a hand closing around a ball. This is called induced fit. Inside the closed pocket, the enzyme’s side chains push and pull on the substrate’s atoms, holding it in a strained shape that is close to the halfway point of the reaction.

The figure shows both ideas, a close-up of a real active site from syn3A’s glycolysis (triose-phosphate isomerase, gene tpiA), and the cycle every enzyme goes through: bind, react, release.

The record holder is an enzyme called orotidine 5′-phosphate decarboxylase. It removes a carbon dioxide from a molecule on the way to making the RNA letter U. Without the enzyme, in water at room temperature, half the molecules would react in about 78 million years. The enzyme speeds the reaction up about 10¹⁷ times, a hundred thousand million million. Divide 78 million years by 10¹⁷ and you get about 25 milliseconds. You will do that division yourself in the Math toolkit below.

Math toolkit: energy barriers and rates

Picture a reaction as a ball that has to roll over a hill. The starting molecules sit in a valley on the left, the products in a lower valley on the right. To get from one to the other, the molecule has to pass through a strained, high-energy shape at the top of the hill, the transition state. The height of the hill above the starting valley is the energy barrier, measured in kilojoules per mole (kJ/mol).

What gets molecules over the top? Heat. Molecules are always being bumped by their neighbors, and now and then a lucky bump gives one enough energy. In Lesson 2 you met the Boltzmann factor: the share of molecules that have at least an energy $E$ is about

$$\text{share} = e^{-E/RT}$$

where $RT$ is the energy of ordinary heat jiggling, per mole. (It is the same as $kT$ per molecule; Lesson 2 used them interchangeably.) At body temperature, 37 °C, $RT = 2.58$ kJ/mol. The reaction’s speed follows that share: twice as many molecules with enough energy means twice as fast.

Now lower the barrier by some amount $x$. The share goes from $e^{-E/RT}$ to $e^{-(E - x)/RT}$, which is $e^{x/RT}$ times bigger. When is that 10 times bigger? When $x/RT = \ln 10 = 2.303$, because $e^{2.303} = 10$. So

$$x = RT \times \ln 10 = 2.58 \times 2.303 = 5.94 \text{ kJ/mol}$$

Three examples, from small to huge:

• Lower the barrier by 5.9 kJ/mol: $10^1$, so 10 times faster.
• Lower it by $3 \times 5.94 = 17.8$ kJ/mol: $10^3$, so 1,000 times faster. Each 5.94 multiplies by another 10, so three of them make $10 \times 10 \times 10$.
• Lower it by $17 \times 5.94 = 101$ kJ/mol: $10^{17}$ times faster. That is the record enzyme’s speed-up. (The record was measured at 25 °C, where $RT$ is a bit smaller, 2.48 kJ/mol, so there the step is 5.71 kJ/mol and the drop is $17 \times 5.71 = 97$ kJ/mol.)

Going the other way: a speed-up of $10^n$ needs a drop of $n \times 5.94$ kJ/mol, and a drop of $x$ gives a speed-up of $10^{x/5.94}$.

Now turn a speed-up into time. In Lessons 3 and 4 you met the rate constant $k$, the chance per second that one molecule reacts, and its link to the half-time:

$$t_{1/2} = \frac{\ln 2}{k} \qquad\text{and}\qquad k = \frac{\ln 2}{t_{1/2}}$$

For the record reaction in water, $t_{1/2}$ = 78 million years. A year has about $3.156 \times 10^7$ seconds, so that is $78 \times 10^6 \times 3.156 \times 10^7 = 2.46 \times 10^{15}$ seconds, and

$$k = \frac{0.693}{2.46 \times 10^{15}\ \text{s}} = 2.8 \times 10^{-16} \text{ per second}$$

The enzyme multiplies $k$ by $10^{17}$: $2.8 \times 10^{-16} \times 10^{17} = 28$ per second. The new half-time is $0.693 \div 28 \approx 0.025$ seconds, the 25 milliseconds from before.

Try it below. The ball is a molecule rattling in the left valley. Turn on the enzyme to open a lower path through the hill, and load the record preset to see the numbers above in the table.

Over the hill

Sound on

A reaction as a hill to climb: the molecule starts in the left valley and must reach the top (the transition state) before it can roll down to the products. Turn the enzyme on, and watch the barrier and the waiting time.

Enzyme offEnzyme on

Show the way back

Orotidine decarboxylaseA peptide bondCustom

To react, a molecule must borrow enough energy from random bumps to get over the top. The higher the hill, the rarer that is.

One molecule waits about 110 million years between crossings. The ball is a cartoon, not to scale.

Barrier without the enzyme161.8 kJ/molThe enzyme lowers it by97 kJ/molTemperature25 °C

kT at 25 °C: RT = 2.48 kJ/molevery 5.71 kJ/mol lower = 10× faster

	No enzyme	With the enzyme
Barrier	161.8 kJ/mol	64.8 kJ/mol
Share with enough energy, e−barrier/RT	4.5 × 10⁻²⁹	4.4 × 10⁻¹²
Rate constant k (per second; rough, Eyring)	2.8 × 10⁻¹⁶	28
Half-time, ln 2 ÷ k (rough)	78 million years	25 ms

97 kJ/mol lower → 10^17.0 times fasterspeed-up e^ΔΔG/RT = 10^x, x = 97 ÷ 5.71 = 17.0

Orotidine decarboxylase. Orotic acid loses CO₂ (a step in making the bases C and U) with a half-time of 78 million years in water at room temperature. The enzyme orotidine 5′-phosphate decarboxylase does it 10¹⁷ times faster: in about 25 milliseconds. Radzicka & Wolfenden 1995; the barriers are worked out from the half-times with Eyring’s formula.

Rates from Eyring’s formula k = (kT/h)·e^{−barrier/RT}, a rough approximation; kT/h = 6.21 × 10¹² per second at 25 °C. RT = 8.314 J/(mol·K) × temperature in kelvin. The ball’s crossings are a cartoon squeezed onto a screen timescale: real waits are in the readouts.

There is one thing an enzyme can’t do. It lowers the barrier for the way back by exactly the same amount, because both directions pass through the same transition state at the top. So if it makes the forward reaction 1,000 times faster, it makes the reverse reaction 1,000 times faster too. The balance point, where forward and backward run equally fast, doesn’t move. An enzyme can’t make a reaction go further, only get to its end sooner. Turn on Show the way back in the widget: $K$, the ratio that sets the balance point, is the same with and without the enzyme.

Check yourselfAt 37 °C, an enzyme lowers a barrier by 23.8 kJ/mol. How much faster does the reaction go?

$23.8 \div 5.94 = 4.0$, so $10^4$: about 10,000 times faster.

Check yourselfA reaction has a half-time of 10 hours. An enzyme speeds it up a million times. What is the new half-time?

10 hours is $10 \times 3{,}600 = 36{,}000$ seconds. Divide by $10^6$: 0.036 seconds, or 36 milliseconds.

Check yourselfAn enzyme speeds up the forward reaction 1,000 times. What does it do to the reverse reaction?

It speeds that up 1,000 times too, because it lowers both barriers by the same amount. Where the reaction ends up doesn’t change.

How fast an enzyme can work

An enzyme has one active site (some have a few) and handles one molecule at a time. Each cycle has three steps: the substrate S binds to the enzyme E, the bound complex ES reacts, and the product P is released:

$$\text{E} + \text{S} \rightleftharpoons \text{ES} \rightarrow \text{E} + \text{P}$$

The double arrow means the substrate can also fall off again before reacting. After the last step the enzyme is free and ready for the next substrate.

So an enzyme has a top speed. Even if a new substrate is waiting the instant the last product leaves, the reaction step itself takes time. The number of molecules one enzyme can finish per second, when it never has to wait, is its turnover number, $k_\text{cat}$ (“k cat”, short for catalytic). Phosphofructokinase, one of syn3A’s glycolysis enzymes, has a $k_\text{cat}$ of about 111 per second (a value measured for E. coli’s version of the enzyme).

The whole cell’s top speed for that reaction, $V_\text{max}$, is $k_\text{cat}$ times the number of enzyme copies. syn3A has about 458 copies of phosphofructokinase, so

$$V_\text{max} = k_\text{cat} \times \text{copies} = 111 \times 458 \approx 50{,}800 \text{ reactions per second}$$

How close an enzyme gets to its top speed depends on how much substrate is around. When substrate is scarce, the enzyme spends most of its time waiting for one to bump into it, and more substrate means more bumps and a faster rate. When substrate is plentiful, the enzyme is always busy: a new molecule arrives the moment the last one leaves. Adding even more then does nothing, because the enzyme is already working flat out. This leveling off is called saturation.

In Enzyme Rush, you are the enzyme. Catch a sugar, press Catalyze, and listen: while you are busy, new sugars bounce off with a thud. Each round puts one point on the chart. Watch the points climb and then flatten out.

3D game · Enzyme Rush

Be the enzyme

Stars 0/18

Sound on

1. Catch and convertCatch2. Turn up the sugarSugar3. Find KmKm4. Look-alikesLook-alikes5. The clampClamp6. Real speedSpeed

[S] 2 mM

Level 1 · Catch and convertSlowed down ~500×

Rate against [S]

your rounds

Start

1. Catch and convert

You are an enzyme: a protein that speeds up one chemical reaction. Orange sugars (glucose 6-phosphate, the first stop of glycolysis) drift around you.

When one lands in the groove between your two halves (the active site), press Catalyze, press Space, or tap the enzyme. Your halves close, the sugar is reshaped into fructose 6-phosphate (lime), and it pops out.

You hold one sugar at a time. While you are busy, others bump off: listen for the thud.

Goal: 30 seconds at [S] = 2 mM. Make 10 products for one star, 13 for two, 16 for three.

The game is slowed down so you can play it: a real enzyme like this one (syn3A’s glucose-6-phosphate isomerase, pgi) does this hundreds of times a second.

Math toolkit: saturation curves

Let’s build the formula for the rate from the game, one piece at a time. Each product takes two stretches of time: the wait until a substrate docks, and the work of reacting it.

• The work time is $1/k_\text{cat}$. In Enzyme Rush, $k_\text{cat}$ = 2 per second, so each reaction takes 0.5 seconds.
• The wait time depends on how much substrate there is. Substrate docks at a rate $k_\text{on} \times [\text{S}]$, where [S] is the substrate concentration and $k_\text{on}$ is a binding rate constant, as in Lesson 4. The average wait is one over that. In the game, $k_\text{on}$ = 0.4 per second per mM, so at 5 mM substrate molecules dock $0.4 \times 5 = 2$ times a second, and the wait is 0.5 seconds.

At 5 mM in the game, one product takes $0.5 + 0.5 = 1$ second: a rate of 1 per second, exactly half the top speed of 2. Try two more:

[S]	Wait	Work	Time per product	Rate	Share of top speed
0.56 mM	4.5 s	0.5 s	5 s	0.2 per s	10%
5 mM	0.5 s	0.5 s	1 s	1 per s	50%
45 mM	0.056 s	0.5 s	0.556 s	1.8 per s	90%

The wait shrinks as you add substrate, but the work never does. That fixed half second is why the rate can never pass 2 per second.

Now the same thing with letters. The rate $v$ of one enzyme is one over the time per product:

$$v = \frac{1}{\text{wait} + \text{work}} = \frac{1}{\dfrac{1}{k_\text{on}[\text{S}]} + \dfrac{1}{k_\text{cat}}}$$

Fractions inside fractions are hard to read, so multiply the top and the bottom by $k_\text{on}[\text{S}] \times k_\text{cat}$. That doesn’t change the value:

$$v = \frac{k_\text{cat}\, k_\text{on} [\text{S}]}{k_\text{cat} + k_\text{on}[\text{S}]}$$

Then divide the top and the bottom by $k_\text{on}$:

$$v = \frac{k_\text{cat} [\text{S}]}{k_\text{cat}/k_\text{on} + [\text{S}]}$$

The group $k_\text{cat}/k_\text{on}$ is a concentration. It gets its own name, the Michaelis constant, $K_\text{m}$. In the game it is $2 \div 0.4 = 5$ mM. Multiply by the number of enzymes to get the whole cell’s rate, and you have the Michaelis–Menten equation:

$$v = V_\text{max} \, \frac{[\text{S}]}{K_\text{m} + [\text{S}]}$$

Read it as two parts. $V_\text{max}$ is the top speed, and the fraction $[\text{S}]/(K_\text{m} + [\text{S}])$ is the share of the time each enzyme is busy. (A real enzyme can also let go of its substrate without reacting. That changes what $K_\text{m}$ is made of, but not the shape of the curve.)

Every Michaelis–Menten curve has the same shape, whatever the enzyme. Measure [S] in units of $K_\text{m}$: if $[\text{S}] = n \times K_\text{m}$, the share of top speed is $n/(n + 1)$.

[S]	Share of top speed
$K_\text{m}/9$	$\tfrac{1/9}{1/9 + 1} = \tfrac{1}{10}$ = 10%
$K_\text{m}$	$\tfrac{1}{2}$ = 50%
$3 K_\text{m}$	$\tfrac{3}{4}$ = 75%
$9 K_\text{m}$	$\tfrac{9}{10}$ = 90%

Going from half speed to 90% takes nine times as much substrate. Getting the last 10% takes forever.

You have seen this curve before. In Lesson 4, the fraction of a protein bound to its partner was $[\text{L}]/(K_\text{d} + [\text{L}])$: half bound when the partner’s concentration equals $K_\text{d}$. The Michaelis–Menten fraction is the same shape, with $K_\text{m}$ in place of $K_\text{d}$. That makes sense: an enzyme is busy exactly when it has a substrate bound.

Now write it in Python. The function takes the substrate concentration, $V_\text{max}$ and $K_\text{m}$, and returns the rate. One check uses syn3A’s real glucose-6-phosphate isomerase.

Python exercise

How fast is an enzyme? Michaelis–Menten

from virtualcell.metabolic_model import KINETICS def michaelis_menten(s, vmax, km): """The rate of an enzyme at substrate concentration s (the Michaelis–Menten equation).""" # s: how much substrate there is, in mM # vmax: the top speed, reached when every enzyme molecule is busy all the time # km: the substrate concentration that gives exactly half the top speed, in mM # # Step 1: the share of the time each enzyme is busy is s / (km + s). # Check it in your head: s = km gives km / (km + km) = 1/2. # Step 2: the rate is vmax times that share. # Return the rate. Careful with brackets: Python does × and ÷ before +. raise NotImplementedError("Write your code here") # Remove the #s below to try it. syn3A's glucose-6-phosphate isomerase (PGI) # has Km = 0.28 mM and a top speed (kcat) of 929 reactions per second. # pgi = KINETICS["PGI"] # for s in [0.03, 0.28, 1, 2.5, 10]: # print(f"{s:5} mM G6P: {michaelis_menten(s, pgi['kcat'], pgi['km']):6.1f} per second")

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfAn enzyme has Km = 0.2 mM and Vmax = 100 per second. How fast does it go at 0.2 mM? At 1.8 mM?

At 0.2 mM, [S] = $K_\text{m}$, so half speed: 50 per second. At 1.8 mM, [S] = $9 K_\text{m}$, so 90%: 90 per second.

Check yourselfHow much substrate does an enzyme need to run at 75% of its top speed?

$3 K_\text{m}$. In general, for a share $f$ you need $[\text{S}] = K_\text{m} \times f/(1 - f)$, and $0.75/0.25 = 3$.

Check yourselfWhen [S] is much smaller than Km, what happens to the rate if you double [S]?

It almost doubles. With [S] tiny, $K_\text{m} + [\text{S}] \approx K_\text{m}$, so $v \approx V_\text{max} [\text{S}]/K_\text{m}$, a straight line through zero. The enzyme is mostly waiting, and twice the substrate means twice the bumps.

Real enzymes, real numbers

What do real enzymes look like? Arren Bar-Even, Ron Milo and colleagues collected measurements of several thousand enzymes in 2011. The typical (median) enzyme has a $k_\text{cat}$ of about 10 per second and a $K_\text{m}$ of about 0.1 mM. Enzymes of central metabolism, such as those of glycolysis, are faster, about 30 per second. The champions are carbonic anhydrase and superoxide dismutase, at a million to ten million reactions per second; restriction enzymes, which cut DNA, manage only one every 10 to 100 seconds.

Lesson 4 showed that two molecules can’t meet faster than diffusion brings them together: about $10^9$ per second per molar. The ratio $k_\text{cat}/K_\text{m}$ says how well an enzyme does when its substrate is scarce, and it has the same units. For the typical enzyme it is $10 \div 0.0001\ \text{M} = 10^5$ per second per molar, about 10,000 times below that limit. A few enzymes, triose-phosphate isomerase among them, come close to the limit: almost every substrate that touches them reacts.

Here are five of syn3A’s own glycolysis enzymes. The $K_\text{m}$ values are the ones syn3A’s 2022 whole-cell model uses, and the $k_\text{cat}$ values are the ones Breuer and colleagues used. Both were measured in related enzymes, used in syn3A’s model: nobody has measured them in syn3A itself. The copy numbers come from counting syn3A’s proteins.

Enzyme (gene)	Substrate	$K_\text{m}$ (mM)	$k_\text{cat}$ (per s)	Copies per cell	Top speed, $k_\text{cat}$ × copies (per s)
Glucose-6-phosphate isomerase (pgi)	G6P	0.28	929	266	247,000
Phosphofructokinase (pfkA)	F6P	0.21	111	458	50,800
Pyruvate kinase (pyk)	PEP	0.075	3,204	551	1,765,000
Lactate dehydrogenase (ldh)	pyruvate	0.96	643	1,100	707,000
Aldolase (fbaA)	FBP	0.005	2.95	227	670

In Km Detective you work out $K_\text{m}$ and $V_\text{max}$ the way a biochemist does: from measurements. Each case gives you dots, one per experiment, and you fit a curve over them. The data are simulated, with noise, from the real values in the table.

Game · Km Detective

Case 1. Read the curve

Stars 0/18

Sound on

1. Read the curve2. The first Km, 19133. Measure in the right place4. A noisy bench5. The fructose case6. The line-up

Case 1

A made-up enzyme

Substrate: substrate · Rate: µM per second

Lay the curve over the dots: set the top speed (Vmax), then Km.

Made-up numbers, with no noise: the dots sit exactly on the true curve.

Where the numbers come from

A made-up enzyme for practice, so the numbers are round.

Measured your curve

Case 1: Read the curve

Each dot is one experiment: mix the enzyme with some substrate and measure how fast product appears.

More substrate, faster reaction, up to a ceiling: the top speed, Vmax. Km is the substrate concentration that gives half the top speed.

Drag the gold handle up or down to set Vmax, and the indigo handle left or right to set Km. The sliders under the chart do the same.

Open the case

Vmax

6 µM/s

−+

Km

9 mM

−+

log scale: each step multiplies

Fit meter 0%★★★★★★

How close your curve is to the best least-squares fit (the curve with the smallest total squared miss).

x-axis: even

Simulated data from published constants. syn3A enzyme values: Metabolic model iMB155 from Breuer et al. 2019, eLife 8:e36842, CC BY 4.0. Km values from Thornburg et al. 2022 (Cell 185:345). Invertase: Michaelis & Menten 1913, reanalyzed by Johnson & Goody 2011.

Now compare an enzyme’s top speed with the work the cell needs from it. syn3A takes up about 15,000 glucose molecules per second, and every one passes through phosphofructokinase. Its 458 copies can do about 50,800 per second, so

$$\text{busy share} = \frac{\text{flow needed}}{\text{top speed}} = \frac{15{,}000}{50{,}800} \approx 0.30$$

Phosphofructokinase works at about 30% of its top speed. That is typical: in Breuer’s model, most enzymes could go 10 to 1,000 times faster than the flow the cell needs. The spare capacity lets the cell speed up when it has to.

Now aldolase. Its borrowed $k_\text{cat}$ of 2.95 per second, times 227 copies, gives a top speed of only 670 per second. But aldolase also has to handle all 15,000 sugars a second. That is about 20 times more than it could do. The cell is alive, so the borrowed number must be wrong for syn3A: perhaps measured in a different organism or under different conditions. Checking capacity against flow is a quick way to catch such errors.

Two short functions: an enzyme’s top speed from $k_\text{cat}$ and its copies, and how busy it is for a given flow. The success message prints all five enzymes against the flow syn3A’s model needs.

Python exercise

An enzyme's top speed in the cell

from virtualcell.constants import SYN3A_GLUCOSE_PER_S # 15,000 glucose taken up per second from virtualcell.metabolic_model import KINETICS def max_rate(kcat, copies): """The top speed of all the copies of an enzyme in one cell, in reactions per second.""" # kcat: reactions per second that ONE enzyme molecule manages at full speed # copies: how many molecules of that enzyme the cell has # All the copies work side by side, so their top speeds add up. raise NotImplementedError("Write your code here") def busy_fraction(flux_per_second, kcat, copies): """How busy each copy must be to carry flux_per_second, as a fraction (0.3 means 30%).""" # Step 1: the top speed of all the copies together is max_rate(kcat, copies). # Step 2: the share of that top speed the cell uses is the flux it needs # divided by the top speed. # Return a fraction (such as 0.3), not a percent (30). raise NotImplementedError("Write your code here") # Remove the #s below to try it on syn3A's phosphofructokinase (PFK), which # handles every glucose the cell takes up. # pfk = KINETICS["PFK"] # print(f"{max_rate(pfk['kcat'], pfk['copies']):,} reactions per second at most") # print(f"busy {busy_fraction(SYN3A_GLUCOSE_PER_S, pfk['kcat'], pfk['copies']):.1%} of the time")

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Inhibitors: slowing enzymes down

An inhibitor is a molecule that slows an enzyme down. Cells use them to control their pathways, and many drugs are inhibitors. They come in two main kinds, and you can tell them apart by what they do to the saturation curve.

A competitive inhibitor looks like the substrate. It fits into the active site but can’t react, so while it sits there, the enzyme is blocked. Substrate and inhibitor compete for the same spot. Add enough substrate and it wins most of the time, so the top speed $V_\text{max}$ doesn’t change. But the enzyme needs more substrate to reach half speed: the apparent $K_\text{m}$ goes up.

$$K_\text{m}^\text{apparent} = K_\text{m} \times \left(1 + \frac{[\text{I}]}{K_\text{i}}\right)$$

[I] is the inhibitor’s concentration, and $K_\text{i}$ is how tightly it binds: the inhibitor concentration that occupies the enzyme half the time, like Lesson 4’s $K_\text{d}$. At $[\text{I}] = K_\text{i}$, the bracket is 2, so $K_\text{m}$ doubles.

A noncompetitive inhibitor binds somewhere else on the enzyme and locks it, whether or not a substrate is in the pocket. More substrate doesn’t help, because the substrate isn’t what’s being pushed out. A locked enzyme is simply out of action. So $V_\text{max}$ goes down, and $K_\text{m}$ stays the same:

$$V_\text{max}^\text{apparent} = \frac{V_\text{max}}{1 + [\text{I}]/K_\text{i}}$$

At $[\text{I}] = K_\text{i}$, half the enzymes are locked and the top speed halves.

Leonor Michaelis and Maud Menten wrote down the saturation curve in 1913, studying invertase, an enzyme that splits table sugar (sucrose) into glucose and fructose. They followed the reaction by shining polarized light through the sugar solution, because sucrose and its products twist the light by different amounts. Re-analyzed in 2011 by Kenneth Johnson and Roger Goody, their data give $K_\text{m}$ = 16.7 mM and $V_\text{max}$ = 0.76 mM per minute.

They also saw that the enzyme’s own products slow it down. Fructose is a competitive inhibitor with $K_\text{i}$ = 58.8 mM, and glucose one with $K_\text{i}$ = 91 mM. With 100 mM fructose in the mix:

$$K_\text{m}^\text{apparent} = 16.7 \times \left(1 + \frac{100}{58.8}\right) = 16.7 \times 2.70 = 45.1 \text{ mM}$$

So with 16.7 mM sucrose, the enzyme drops from half speed (0.38 mM per minute) to $0.76 \times 16.7/(45.1 + 16.7) = 0.21$ mM per minute.

Cells use inhibition to control their pathways. Often the final product of a pathway inhibits one of its first enzymes, so the pathway slows down by itself when enough product has piled up. This is called feedback inhibition, and Lesson 6 builds whole circuits out of feedback.

Write both kinds of inhibited rate. The checks use invertase’s real constants.

Python exercise

Inhibitors: two ways to slow an enzyme

def competitive_rate(s, vmax, km, inhibitor, ki): """The rate with a competitive inhibitor, which blocks the active site.""" # A competitive inhibitor looks like the substrate and sits in the active # site. It makes the enzyme seem to need more substrate: the apparent Km is # km × (1 + inhibitor / ki) # where inhibitor is its concentration (mM) and ki says how tightly it binds # (mM: the smaller, the tighter). vmax doesn't change. # Step 1: work out the apparent Km. # Step 2: return the Michaelis–Menten rate with vmax and the apparent Km. # (michaelis_menten, below, is the same as your answer to exercise 1.) raise NotImplementedError("Write your code here") def noncompetitive_rate(s, vmax, km, inhibitor, ki): """The rate with a noncompetitive inhibitor, which switches enzymes off from elsewhere.""" # A noncompetitive inhibitor binds away from the active site and stops the # enzyme working. It lowers the top speed: the apparent vmax is # vmax / (1 + inhibitor / ki) # and Km doesn't change. # Step 1: work out the apparent vmax. # Step 2: return the Michaelis–Menten rate with the apparent vmax and km. raise NotImplementedError("Write your code here") def michaelis_menten(s, vmax, km): """The same as your answer to exercise 1.""" return vmax * s / (km + s) # Remove the #s below to try it on Michaelis and Menten's own enzyme, invertase, # which splits sucrose: Km 16.7 mM, top speed 0.76 mM per minute. Fructose, one # of the two sugars it makes, is a competitive inhibitor with Ki = 58.8 mM. # print("sucrose rate (mM per minute)") # print(" (mM) alone + 100 mM fructose") # for sucrose in [5, 16.7, 50, 200, 1000]: # plain = michaelis_menten(sucrose, 0.76, 16.7) # slowed = competitive_rate(sucrose, 0.76, 16.7, 100, 58.8) # with 100 mM fructose # print(f"{sucrose:7} {plain:.3f} {slowed:.3f}")

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

This one is a challenge: given a list of measured rates, find the $K_\text{m}$ and $V_\text{max}$ that fit them best, as you did by hand in Km Detective. The data are simulated from real enzymes, with noise added. The starter walks you through a search.

Python exercise

Challenge: fit Km and Vmax to data

def fit_michaelis_menten(s_values, rates): """The (vmax, km) whose Michaelis–Menten curve best fits the measured rates.""" # "Best" means the smallest sum_of_squares (below): the curve that misses # the measured points by the least, with each miss squared so that misses # above and below the curve both count. # You need two numbers. A good plan: # Step 1: for one fixed km, the best vmax has a formula. With # f = s / (km + s) for each point, # the best vmax is sum(rate × f) / sum(f × f), adding over all points. # (best_vmax, below, does this for you.) # Step 2: try many km values spread out on a log scale, each a fixed factor # bigger than the one before, from min(s_values) / 100 up to # max(s_values) * 100. 101 values: low * (high / low) ** (i / 100) for # i in range(101). For each km, work out sum_of_squares(best_vmax(km, ...), km, ...). # Step 3: keep the km with the smallest sum. Then zoom in: make a new, narrower # range around it (for example from km / 1.5 to km × 1.5) and repeat # step 2 a few times. # Return (vmax, km), with vmax first. raise NotImplementedError("Write your code here") def best_vmax(km, s_values, rates): """For a fixed km, the vmax that gives the smallest sum of squares.""" f = [s / (km + s) for s in s_values] return sum(r * x for r, x in zip(rates, f)) / sum(x * x for x in f) def sum_of_squares(vmax, km, s_values, rates): """How badly the curve (vmax, km) misses the measured rates: the misses squared, added up.""" total = 0.0 for s, rate in zip(s_values, rates): miss = rate - michaelis_menten(s, vmax, km) total += miss * miss return total def michaelis_menten(s, vmax, km): """The same as your answer to exercise 1.""" return vmax * s / (km + s) # Made-up data, simulated from the values syn3A's model uses for lactate # dehydrogenase (kcat 643 per second, Km 0.96 mM for pyruvate), with 10% random noise. LDH_PYRUVATE_MM = [0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4] LDH_RATE_PER_S = [52.4, 138.8, 230.1, 275.7, 422.2, 396.9, 637.7] # Remove the #s below to try it. The data were made with 643 and 0.96: how close do you get? # vmax, km = fit_michaelis_menten(LDH_PYRUVATE_MM, LDH_RATE_PER_S) # print(f"vmax = {vmax:.0f} per second, Km = {km:.3f} mM")

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfA competitive inhibitor is present at a concentration equal to its Ki. What happens to Km and Vmax?

$K_\text{m}$ doubles, because $1 + K_\text{i}/K_\text{i} = 2$. $V_\text{max}$ stays the same: with enough substrate, the enzyme still reaches its top speed.

Check yourselfAdding lots of substrate doesn't bring back an inhibited enzyme's top speed. Which kind of inhibitor is it?

Noncompetitive. It locks the enzyme away from the active site, so substrate can’t push it out.

Glycolysis: sugar to energy

A pathway is a chain of reactions where the product of one enzyme is the substrate of the next. Glycolysis (“sugar splitting”) is the pathway that breaks one glucose, a sugar with six carbon atoms, into two molecules of pyruvate, with three carbons each. It has ten steps, in three acts:

1. Invest. The cell spends ATP to stick phosphate groups onto the sugar. That makes it more reactive, and the charged phosphates keep it from leaking out of the cell.
2. Split. The six-carbon sugar is cut into two three-carbon halves.
3. Pay off. Each half is rearranged, step by step, in a way that releases energy, and that energy is stored as ATP.

One of the pay-off steps also strips electrons from the sugar. They are picked up by a carrier molecule, NAD⁺ (nicotinamide adenine dinucleotide), which becomes NADH when loaded. Think of NAD⁺ as an empty cart and NADH as a full one. Keep that cart in mind: it will turn out to be a problem.

Count the textbook version per glucose. The invest steps spend 2 ATP. The pay-off steps make 2 ATP per half, and there are two halves, so 4 ATP. Net:

$$-2 + 4 = +2 \text{ ATP per glucose}$$

plus 2 NADH and 2 pyruvate.

Here is the whole pathway as syn3A runs it, with its real gene names and the length of each protein. Orange badges mark six-carbon sugars, lime ones the three-carbon halves, and violet the two-carbon pieces. Follow the arrows down from glucose to pyruvate, then out to lactate or acetate at the bottom.

Glycolysis in syn3A

The ten steps in syn3A

Here are the steps with the genes that do them in syn3A, from its GenBank record. The locus tags (JCVISYN3A_ and four digits) are the genes’ IDs; you met them in Lesson 3.

Step	What happens	Enzyme	Gene	Locus tag
In	glucose comes in and gets a phosphate: glucose 6-phosphate (G6P)	the PTS (four proteins)	ptsI, ptsH, crr, ptsG	0233, 0694, 0234, 0779
1	G6P is rearranged into fructose 6-phosphate (F6P)	glucose-6-phosphate isomerase	pgi	0445
2	a second phosphate, from ATP: fructose 1,6-bisphosphate (FBP)	phosphofructokinase	pfkA	0220
3	FBP is cut into two three-carbon halves, DHAP and G3P	aldolase	fbaA	0131
4	DHAP is turned into a second G3P	triose-phosphate isomerase	tpiA	0727
5	G3P gets a phosphate and gives two electrons to NAD⁺: 1,3-bisphosphoglycerate	glyceraldehyde-3-phosphate dehydrogenase (GAPDH)	gapDH	0607
6	one phosphate moves to ADP, making ATP: 3-phosphoglycerate (3PG)	phosphoglycerate kinase	pgk	0606
7	the other phosphate moves along the molecule: 2PG	phosphoglycerate mutase	pgm	0729
8	a water comes off: phosphoenolpyruvate (PEP)	enolase	eno	0213
9	the last phosphate moves to ADP, making ATP: pyruvate	pyruvate kinase	pyk	0221

Everything from step 5 on happens twice per glucose, once for each half.

Textbooks count ten steps, starting with hexokinase, an enzyme that puts a phosphate from ATP onto glucose. syn3A has no hexokinase gene: it was removed when the genome was cut down. So the table has the PTS, then nine enzymes.

Instead, glucose enters through the PTS (phosphotransferase system), four proteins that work like a bucket chain. ptsI takes a phosphate from PEP, the molecule that step 9 turns into pyruvate. It hands the phosphate to ptsH, which hands it to crr, which hands it to ptsG. ptsG sits in the membrane, and it attaches the phosphate to glucose as it carries the sugar into the cell. So glucose arrives already as G6P, and syn3A has paid for it with one of its own PEPs, from the end of its own assembly line.

Does syn3A still make a profit? Count again, per glucose:

• Phosphofructokinase spends 1 ATP. There is no hexokinase to spend another.
• Phosphoglycerate kinase makes 1 ATP per half: 2 ATP.
• Of the two PEPs, one paid for the glucose. Only the other reaches pyruvate kinase: 1 ATP.

$$-1 + 2 + 1 = +2 \text{ ATP per glucose}$$

The same profit as the textbook version, paid in a different way. And both PEPs still end up as pyruvate: the one that paid the PTS becomes pyruvate when it gives up its phosphate.

In Glycolysis Line you build the pathway as a factory line. Put each enzyme in its slot, follow the carbons through the split, watch the ATP wallet, then switch to syn3A’s genes and its PTS.

3D game · Glycolysis Line

Run the sugar factory

Mistakes 0Stars 0/18

Sound on

1. Build the line2. Follow the carbons3. Invest 2, earn 44. syn3A’s way in5. The NAD⁺ crunch6. The acetate bonus

1. Build the line

Glycolysis is a factory line of ten enzymes. Each one makes one small change to the sugar and hands it on.

Each empty slot shows what goes in and what comes out. Tap an enzyme card, then the slot where it belongs. When all ten are in, a glucose runs down the line.

Start building

PlatesTextbooksyn3A

The ten stations · tap one for details

1?→glucose → G6P2?→G6P → F6P3?→F6P → FBP4?→FBP → DHAP + G3P5?→DHAP → G3P6?→G3P → 1,3-BPG7?→1,3-BPG → 3PG8?→3PG → 2PG9?→2PG → PEP10?→PEP → pyruvate

Genes and locus tags from syn3A’s GenBank record; reactions from the metabolic model iMB155 (Breuer et al. 2019, eLife 8:e36842, CC BY 4.0). The molecules are simplified: hydrogens, water and H⁺ are left out, and the line runs far slower than a cell, which makes about 15,000 glucose into pyruvate every second.

Check yourselfHow many ATP does syn3A make from 5 glucose sent all the way to lactate?

2 per glucose, so 10 ATP. It also fills 10 NAD⁺ carts with electrons, and, as the next section explains, the lactate step empties all 10 again.

Check yourselfSix carbons go in with each glucose. How many carbons come out as pyruvate?

Six: two pyruvates of three carbons each. Glycolysis rearranges the carbons but loses none.

Paying back NADH: lactate and acetate

Step 5, GAPDH, needs an empty NAD⁺ cart for every G3P. syn3A’s whole-cell model starts with about 2.18 mM NAD⁺, which in a cell as small as syn3A is about 44,000 molecules. (You will see how to turn mM into molecules in the next Math toolkit.) Glycolysis runs at about 15,000 glucose per second, two G3P each, so GAPDH fills about 30,000 carts a second. If nothing emptied them:

$$\frac{44{,}000 \text{ carts}}{30{,}000 \text{ per second}} \approx 1.5 \text{ seconds}$$

In about a second and a half, every cart would be full and glycolysis would stop.

The fix is lactate dehydrogenase (gene ldh, JCVISYN3A_0475). It hands the electrons back to pyruvate:

$$\text{pyruvate} + \text{NADH} \rightarrow \text{lactate} + \text{NAD}^+$$

Each pyruvate empties exactly one cart, and each glucose filled two carts and made two pyruvates, so the books balance. The lactate leaves the cell. No ATP is made in this step; it exists only to return the carts. Your muscles do the same thing when they run short of oxygen.

syn3A has a second route, which pays better. Pyruvate can be cut into a two-carbon acetyl group plus carbon dioxide (genes pdhC and lpdA, 0227 and 0228), passed to phosphate (pta, 0229), and finally to ADP by acetate kinase (ackA, 0230), which makes one extra ATP per pyruvate. The acetate leaves the cell. But the first step of this route fills another NAD⁺ cart instead of emptying one.

Those carts are emptied by an NADH oxidase, which hands the electrons to oxygen and makes water. Oddly, no gene for it is known in syn3A: the one it had was removed when the genome was cut down. Both syn3A models need the reaction, so they assume some other protein does the job. And acetate can only leave the cell so fast. The model borrows a limit from a related bacterium: at most about 6.9 acetate for every 7.4 glucose. So syn3A sends some pyruvate each way. With all acetate, each glucose could make up to 4 ATP; with the limit, the model’s best is about 2.9.

ATP: the energy currency

All that ATP: what is it, and why does every cell use it? ATP (adenosine triphosphate) is a nucleotide, like the letters of RNA. It has three parts: a base called adenine, a sugar called ribose, and a tail of three phosphate groups. The phosphates are all negatively charged and repel each other. When water breaks the bond to the last one, ATP becomes ADP (adenosine diphosphate) plus a loose phosphate, and energy is released.

How much? Under the standard conditions chemists use for tables, about 30 kJ/mol. But a cell keeps far more ATP than ADP, and that pushes the reaction harder. In E. coli growing on glucose it releases about 47 kJ/mol. As in Lesson 2, we use about 50 kJ/mol, or about 20 kT per ATP: $50 \div 2.58 \approx 19$. One ATP carries

$$\frac{50{,}000 \text{ J/mol}}{6.022 \times 10^{23} \text{ per mol}} = 8.3 \times 10^{-20} \text{ J}$$

The figure shows the molecule, with the bond that breaks; the cycle of ATP being spent and remade; and how small syn3A’s supply is compared with how fast it flows.

Why is it so useful? Because energies add up. Suppose a cell needs a reaction that is uphill: say it needs 20 kJ/mol put in, so on its own it barely happens. If an enzyme ties it to breaking one ATP, the two together release about $50 - 20 = 30$ kJ/mol, and the combined reaction is downhill. This is called coupling. Phosphofructokinase does it: it takes a phosphate from ATP and puts it on the sugar in one move. So do ribosomes, pumps and almost every builder in the cell. ATP is the currency; glycolysis earns it, and everything else spends it.

Math toolkit: pools and flows

A cell’s ATP is a pool: an amount that is constantly drained and refilled. Think of a bathtub with the tap on and the plug out. How long does one molecule stay in the tub, on average? Divide the amount by the flow:

$$\text{turnover time} = \frac{\text{pool}}{\text{flow}}$$

A tub holding 100 liters, with 10 liters a minute flowing through, replaces its water every 10 minutes. It is Lesson 3’s steady state again, read the other way: there, the number of copies was the production rate times the average lifetime.

To use it, we need pools as numbers of molecules. Concentrations come in mM, so here is the conversion for syn3A, one step at a time.

• 1 mM means $10^{-3}$ moles per liter.
• syn3A is a sphere 200 nm in radius, with a volume of 0.0335 femtoliters (fL), and 1 fL is $10^{-15}$ liters. So the volume is $3.35 \times 10^{-17}$ liters.
• Moles in the cell: $10^{-3} \times 3.35 \times 10^{-17} = 3.35 \times 10^{-20}$.
• Molecules: multiply by Avogadro’s number from Lesson 1: $3.35 \times 10^{-20} \times 6.022 \times 10^{23} \approx 20{,}170$.

Now the ATP pool. The whole-cell model starts syn3A with 3.65 mM ATP, a value estimated from other bacteria (nobody has measured it in syn3A). That is $3.65 \times 20{,}170 \approx$ 74,000 ATP molecules. The same model makes about 35,000 ATP per second, and spends about as many. So

$$\text{turnover time} = \frac{74{,}000}{35{,}000 \text{ per second}} \approx 2.1 \text{ seconds}$$

The cell’s entire supply of ATP is used up and remade about every 2 seconds. Other pools are even faster. PEP, at 0.041 mM, is about 830 molecules, made at about 30,000 per second: replaced every 30 milliseconds or so.

Metabolic models like iMB155 don’t count molecules per second. They measure a flux, the speed of a reaction, in millimoles per gram of dry weight per hour (mmol per gram dry weight per hour), which lets the same numbers work for a test tube of cells of any size. To turn that into molecules per second in one syn3A cell:

• 1 mmol is $10^{-3}$ moles.
• One syn3A cell has about 10.2 femtograms of dry mass, $10.2 \times 10^{-15}$ grams (a 400 nm sphere, mostly water, the rest protein, RNA, DNA and lipid).
• So per cell, per hour: $10^{-3} \times 10.2 \times 10^{-15} = 1.02 \times 10^{-17}$ moles.
• Molecules per hour: $1.02 \times 10^{-17} \times 6.022 \times 10^{23} = 6.14 \times 10^6$.
• Per second: divide by 3,600 seconds in an hour: $6.14 \times 10^6 \div 3{,}600 \approx 1{,}706$.

Check yourselfA pool of 500 molecules is made and used at 1,000 per second. What is its turnover time?

$500 \div 1{,}000 = 0.5$ seconds.

Check yourselfThe model says acetate leaves at 6.9 mmol per gram dry weight per hour. How many acetate molecules is that per second, in one syn3A cell?

$6.9 \times 1{,}706 \approx 11{,}800$ per second.

Check yourselfsyn3A's model starts with about 0.22 mM ADP. About how many molecules is that?

$0.22 \times 20{,}170 \approx 4{,}400$ ADP molecules, about 17 times fewer than ATP.

Where the ATP goes

What does syn3A spend 35,000 ATP a second on? The biggest item is making protein. Each amino acid costs about 4 ATP (Lesson 3), and about 500 ribosomes add about 5 amino acids a second each:

$$500 \times 5 \times 4 = 10{,}000 \text{ ATP per second}$$

Next is importing food: the whole-cell model assumes 1 ATP per molecule brought in, about 5,000 a second. Copying and reading genes takes about 2,900 more. Upkeep takes about 1,000: the membrane leaks protons, and the ATP synthase runs backwards to pump them out again, using at least 0.575 mmol per gram dry weight per hour of ATP, or $0.575 \times 1{,}706 \approx 1{,}000$ a second. Everything else the cell builds, such as lipids and nucleotides, uses the rest.

The published model splits the bill its own way, and one surprise is how much it can’t explain. Of all the ATP it makes, 31% goes to a lump of growth costs that the modelers could not put a name to, 30% to making protein, and 26% to phosphofructokinase, glycolysis’s own investment step. Nucleotides, transport, RNA and proton pumping take 2–3% each.

The ATP Economy runs syn3A’s budget live. The gold tank is the ATP pool. Glycolysis pours coins in; ribosomes, transporters, polymerases, the proton pump and a building crane take them out. Move the sliders and see how fast the tank responds.

Simulation · The ATP Economy

Earn it, spend it, every second

Stars 0/15

Sound on

1. Break even2. Two seconds3. E. coli speed4. The death spiral5. Where does it go?

Challenge: Set the glucose uptake, press Play, and hold ATP between 3 and 4 mM for 20 seconds. Use as little sugar as you can.

74,000 ATP3.67 mM

turnover 0.0 sclock 0.0 s

Where the ATP goes

BudgetFBA view

• Translation37%
• Imports19%
• Gene work11%
• Upkeep4%
• Building30%

1. Break even

This is the ATP budget of one syn3A cell, in real time. Glycolysis on the left turns glucose into ATP (gold coins). Everything else in the cell spends it.

Right now the cell takes up only 6,000 glucose a second, and each glucose pays 2 ATP. The spenders want about 26,900 ATP a second.

Find the smallest uptake that keeps the tank between 3 and 4 mM for 20 seconds. If the tank runs dry, press Reset.

Set up the cell

★★★ 14,500 a second or less · ★★ up to 17,000 · ★ any uptake that holds

A budget model: real cells regulate every one of these flows. The numbers are per second, for one cell, in real time.

PlayReset

×1 real time×¼ slow

Glucose uptake

6,000 a second

−+

2 ATP per glucose: up to 12,000 ATP a second. syn3A: about 15,000.

Pyruvate sent to acetate

0%

−+

Fixed for this challenge. All pyruvate becomes lactate. Each pyruvate sent to acetate pays 1 more ATP (limit 11,800 a second).

Ribosomes

500

−+

Fixed for this challenge. Translation: 500 × 5 amino acids a second × 4 ATP = 10,000 a second. syn3A: about 500.

Ribosome speed

5 amino acids a second

−+

Fixed for this challenge. syn3A: about 5. E. coli: about 10 to 20.

Building

8,000 ATP a second

−+

Fixed for this challenge. Lipids, nucleotides and the costs nobody has itemized (the rest of the budget). Default 8,000.

Imports and fixed costs

Imports on: 5,000 a second

Always on: gene work 2,900 a second (RNA and DNA polymerases) and upkeep 1,000 a second (the ATP synthase pumping protons out). Imports are fixed for this challenge.

The arithmetic, right now

Income with plenty of ATP

2 × 6,000 glucose = 12,000 a second

Spending with plenty of ATP

10,000 + 5,000 + 2,900 + 1,000 + 8,000 = 26,900 a second

Pool in molecules

3.67 mM × 20,170 per mM = 74,000

Turnover time

pool ÷ spending = 74,000 ÷ 0 = 0.0 s

At full speed glycolysis could pay 12,000 a second, but it slows when ATP runs low (it must invest ATP first) or when ADP runs short (ATP is made from ADP). Spenders slow a little as ATP gets scarce. So the pool settles where income equals spending.

Where the numbers come from

• Pool: 1 mM in syn3A = 10⁻³ mol/L × 3.35 × 10⁻¹⁷ L × 6.022 × 10²³ = 20,170 molecules. The whole-cell model starts with 3.65 mM ATP, about 74,000 molecules (Thornburg et al. 2022, Cell 185:345; estimated from E. coli and M. pneumoniae, not measured in syn3A). ATP + ADP is held at 80,000 (about 4 mM).
• Glucose: about 15,000 taken up a second; glycolysis pays 2 ATP each (Thornburg et al. 2022).
• Acetate: at most 6.9 mmol per gram dry weight per hour in syn3A’s model (Breuer et al. 2019, borrowed from M. pneumoniae) × 1,706 ≈ 11,800 a second.
• Translation: 500 ribosomes × 5 amino acids a second × 4 ATP = 10,000 a second (Thornburg et al. 2019; 503 ribosomes counted in a syn3A tomogram).
• Imports: about 5,000 a second, 1 ATP per molecule brought in (Thornburg et al. 2022). Gene work: 18 million ATP per cell cycle ÷ 6,300 s ≈ 2,900 a second (Thornburg et al. 2019, derived).
• Upkeep: at least 0.575 mmol per gram dry weight per hour for proton pumping by the ATP synthase (Breuer et al. 2019) × 1,706 ≈ 1,000 a second. Building: the rest, 8,000.
• Assumed (chosen so the defaults match the model’s pool and flows): glycolysis runs at N² ÷ (6,000² + N²) of full speed and needs ADP (× ADP ÷ (1,000 + ADP)); each spender runs at N ÷ (4,000 + N) of its set rate (Km about 0.2 mM, between the model’s 0.12 for PFK and 0.37 for PGK). Steps of 5 ms.
• FBA view: Metabolic model iMB155 from Breuer et al. 2019, eLife 8:e36842, CC BY 4.0.

How much power is that? If glycolysis makes about 30,000 ATP a second (15,000 glucose to lactate, 2 each), the cell runs on $30{,}000 \times 8.3 \times 10^{-20}\ \text{J} \approx 2.5 \times 10^{-15}$ watts. A 10-watt light bulb uses as much power as about four million billion syn3A cells.

The stoichiometric matrix

Now step back and look at the whole factory at once. Glycolysis alone has a dozen reactions, and syn3A’s model has 338. To reason about all of them together, we write them in one table.

Start with one reaction, phosphofructokinase (PFK):

$$\text{F6P} + \text{ATP} \rightarrow \text{FBP} + \text{ADP}$$

Write it as a column of numbers, one per molecule: how many it uses (negative) and how many it makes (positive). PFK uses one F6P and one ATP, and makes one FBP and one ADP. So its column has −1 for F6P, −1 for ATP, +1 for FBP, +1 for ADP, and 0 for every other molecule. These numbers are the reaction’s stoichiometry (“stoy-key-OM-etry”, the counting of molecules in a reaction).

Now put every reaction’s column side by side. You get a table with one row per metabolite and one column per reaction: the stoichiometric matrix, written $S$. Each entry says what one reaction does to one molecule. Read down a column to see a reaction. Read across a row to see every reaction that makes or uses one molecule.

This lesson uses a core of syn3A’s model: 19 reactions for the PTS, glycolysis, the lactate and acetate routes, the oxidase, and the exits. They are iMB155’s own reactions, with a few pairs of steps joined into one, plus one added reaction, ATPM, that stands for all the cell’s ATP spending. They involve 19 balanced metabolites. A few small molecules, H⁺, water, phosphate, carbon dioxide, oxygen and the glucose outside, are left out of the table, as if the medium supplied or removed them freely.

Here is the matrix. Tap a reaction to read its column, or tap a molecule’s row to see every reaction that touches it. Start with the small toy network, then switch to syn3A’s core.

The reaction ledger

Sound on

Write every reaction as a column of numbers and stack the columns side by side: that table is the stoichiometric matrix, S. Multiply it by the reaction speeds, v, and each row tells you whether that molecule piles up or runs out.

Toysyn3A coreFill it in

Ainv = 4outside→AR1v = 4A→BR2v = 2B→CR3v = 2B→DCoutv = 2C→outsideDoutv = 2D→outside

moleculev =	Ain4	R14	R22	R32	Cout2	Dout2	S·v
A	+1	−1					0
B		+1	−1	−1			0
C			+1		−1		0
D				+1		−1	0

Steady state: every row of S·v is 0. Nothing piles up, nothing runs out.

A ledger for four molecules

Each column is a reaction: negative for what it uses, positive for what it makes, blank for everything else. Each row is a molecule. Tap a reaction or a molecule.

The last column multiplies each row by the speeds v and adds up: how fast that molecule piles up (+) or runs out (−). All zeros is a steady state.

Reaction speeds, v (made-up units)

Ain A comes in4R1 A turns into B4R2 B turns into C2R3 B turns into D2Cout C leaves2Dout D leaves2

The toy network is made up. syn3A core: 19 reactions of glycolysis and its two exits, taken from the metabolic model iMB155 (three are lumped together from several model steps; ATP use is added to stand for everything the cell spends). Metabolic model iMB155 from Breuer et al. 2019, eLife 8:e36842, CC BY 4.0.

Math toolkit: vectors and matrices

This toolkit starts from zero. Two new words: vector and matrix.

A vector is a list of numbers in a fixed order. For example, the speeds of two reactions, PFK at 7.4 and aldolase at 7.4, make the vector $(7.4,\ 7.4)$. A list of every reaction’s speed is called the flux vector, written $v$. The order matters: the first number always belongs to the first reaction, and so on. In Python it is just a list.

The one operation you need is the dot product: take two lists of the same length, multiply the numbers in matching positions, and add the results. For $(2, 3)$ and $(4, 5)$:

$$(2, 3) \cdot (4, 5) = 2 \times 4 + 3 \times 5 = 8 + 15 = 23$$

With three numbers: $(1, -1, 0) \cdot (5, 5, 5) = 5 - 5 + 0 = 0$. That’s all it is: pair up, multiply, add. In Python:

sum(x * y for x, y in zip([2, 3], [4, 5]))   # 23

A matrix is a table of numbers, with rows and columns. Its size is written rows × columns: a 2 × 3 matrix has 2 rows and 3 columns, so 6 entries. Each row of a matrix is itself a list of numbers, a vector.

Matrix times vector means: do a dot product of each row with the vector, one row at a time. Each row gives one number, so the answer is a list with one number per row. Here is a 2 × 2 example:

$$\begin{pmatrix} 1 & -1 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 3 \\ 3 \end{pmatrix} = \begin{pmatrix} 1 \times 3 + (-1) \times 3 \\ 0 \times 3 + 1 \times 3 \end{pmatrix} = \begin{pmatrix} 0 \\ 3 \end{pmatrix}$$

The first row, $(1, -1)$, dotted with $(3, 3)$ gives 0. The second row, $(0, 1)$, gives 3. The vector is written standing up, as a column, but it is the same list.

What does it mean for a metabolic network? Say this matrix is $S$ for two reactions, R1 that makes A from food, and R2 that turns A into B. Row A is $(+1, -1)$: R1 makes A, R2 uses it. Row B is $(0, +1)$: only R2 touches B, and it makes it. With both reactions running at 3 per second, A is made at 3 and used at 3, so it doesn’t change: 0. B is made at 3 and never used, so it piles up: +3 per second. Each entry of $S \cdot v$ is how fast one molecule piles up (positive) or runs out (negative).

Now a 3 × 3. Three molecules in a loop: R1 turns A into B, R2 turns B into C, and R3 turns C back into A.

	R1 (A → B)	R2 (B → C)	R3 (C → A)
A	−1	0	+1
B	+1	−1	0
C	0	+1	−1

With every reaction at 2, $v = (2, 2, 2)$:

• Row A: $(-1)(2) + (0)(2) + (+1)(2) = -2 + 0 + 2 = 0$
• Row B: $(+1)(2) + (-1)(2) + (0)(2) = 0$
• Row C: $(0)(2) + (+1)(2) + (-1)(2) = 0$

Every row is 0: nothing changes. Now speed R1 up to 3, so $v = (3, 2, 2)$. Row A becomes $-3 + 2 = -1$, and row B becomes $3 - 2 = +1$, while C stays 0. A runs down by 1 per second and B piles up by 1 per second.

Now one real row from syn3A’s core. FBP is made by PFK and used by aldolase (the reaction FBA in the model), so its row has +1 under PFK, −1 under FBA, and 0 everywhere else. In the model’s best plan, both run at 7.4 mmol per gram dry weight per hour:

$$(+1)(7.4) + (-1)(7.4) = 0$$

FBP stays level. The whole core matrix has $19 \times 19 = 361$ entries, but only 53 of them are not zero, because each reaction touches only a few molecules. The full model’s matrix is $304 \times 338$, about 100,000 entries, almost all zero.

Build the matrix in Python. Given a list of reactions, each with its stoichiometry as a dict, return the metabolites, the reaction ids, and the rows of $S$. The checks start with a toy and end with syn3A’s core.

Python exercise

Build the stoichiometric matrix

from virtualcell.metabolic_model import METABOLITES, core_reactions def stoichiometric_matrix(reactions): """The stoichiometric matrix of a list of reactions: (metabolites, reaction_ids, S).""" # Each reaction is a dict. reaction["stoich"] says how many of each metabolite # it uses (negative) or makes (positive), for example {"g6p_c": -1, "f6p_c": 1}. # Step 1: metabolites = every metabolite that appears in any reaction's # "stoich", once each, sorted alphabetically. (Collect them in a set, # then use sorted().) # Step 2: reaction_ids = the reactions' "id"s, in the order they were given. # Step 3: S is a list of rows, one row per metabolite. Row i has one number per # reaction: the coefficient of metabolite i in reaction j, or 0 if the # reaction doesn't use it. (stoich.get(name, 0) gives exactly that.) # Return (metabolites, reaction_ids, S). raise NotImplementedError("Write your code here") # A small made-up network to try it on: food A comes in, becomes B, B splits # into C and D, and both of those leave. TOY = [ {"id": "R1", "stoich": {"A": 1}}, # → A {"id": "R2", "stoich": {"A": -1, "B": 1}}, # A → B {"id": "R3", "stoich": {"B": -1, "C": 1, "D": 1}}, # B → C + D {"id": "R4", "stoich": {"C": -1}}, # C → {"id": "R5", "stoich": {"D": -1}}, # D → ] # Remove the #s below to print the toy's matrix (· means 0). # metabolites, ids, S = stoichiometric_matrix(TOY) # print(" " + "".join(f"{r:>4}" for r in ids)) # for name, row in zip(metabolites, S): # print(f"{name:>5}" + "".join(f"{x:>4}" if x else " ·" for x in row)) # print() # # syn3A's core is 19 × 19 and mostly zeros, so print each row as a ledger: # the reactions that make the metabolite (+) and the ones that use it (−). # metabolites, ids, S = stoichiometric_matrix(core_reactions()) # for name, row in zip(metabolites, S): # entries = [f"{rid} {x:+g}" for rid, x in zip(ids, row) if x] # print(f"{METABOLITES[name]['short']:>10}: " + ", ".join(entries))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Then multiply it by a flux vector: one dot product per row. With the model’s best fluxes, all 19 rows of syn3A’s core come out zero.

Python exercise

Mass balance: S times v

from virtualcell.metabolic_model import FBA_CORE_FLUXES, METABOLITES, core_reactions def mass_balance(S, v): """S times v: for each metabolite, how fast it builds up (+) or runs down (−).""" # S is a list of rows, one per metabolite; v is a list of fluxes, one per reaction. # Row i times v: multiply row[j] by v[j] for every j, then add them all up. # That one number says how fast metabolite i builds up (+) or runs down (−). # Do it for every row and return the list of sums, one per metabolite. raise NotImplementedError("Write your code here") def is_steady_state(S, v, tol=1e-6): """True if every metabolite is made exactly as fast as it is used (to within tol).""" # Every entry of mass_balance(S, v) must be close to 0: between -tol and +tol. # (abs(x) <= tol says that. all(...) is True when every item is True.) raise NotImplementedError("Write your code here") def stoichiometric_matrix(reactions): """The same as your answer to exercise 5.""" metabolites = sorted({name for reaction in reactions for name in reaction["stoich"]}) reaction_ids = [reaction["id"] for reaction in reactions] S = [[reaction["stoich"].get(name, 0) for reaction in reactions] for name in metabolites] return metabolites, reaction_ids, S # Remove the #s below to try it. First by hand: 1×3 + (−1)×3 = 0 and 0×3 + 1×3 = 3. # print(mass_balance([[1, -1], [0, 1]], [3, 3])) # should be [0, 3] # # Then syn3A's core with the fluxes that flux balance analysis finds (FBA_CORE_FLUXES): # metabolites, ids, S = stoichiometric_matrix(core_reactions()) # v = [FBA_CORE_FLUXES[r] for r in ids] # for name, change in zip(metabolites, mass_balance(S, v)): # print(f"{METABOLITES[name]['short']:>10}: {round(change, 9) + 0:+.2f}") # + 0 tidies -0.0 into 0.0 # print("steady state:", is_steady_state(S, v))

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfWrite the column for lactate dehydrogenase: pyruvate + NADH → lactate + NAD⁺.

Pyruvate −1, NADH −1, lactate +1, NAD⁺ +1, and 0 for every other molecule.

Check yourselfIn the core's best plan, PGK runs at 14.8, pyruvate kinase at 7.4, acetate kinase at 6.9, PFK at 7.4 and ATPM at 21.7. The first three make ATP and the last two use it. Does ATP balance?

$14.8 + 7.4 + 6.9 - 7.4 - 21.7 = 29.1 - 29.1 = 0$. Yes: ATP is made exactly as fast as it is used.

Check yourselfA matrix has 4 rows and 5 columns. How many entries does it have, and how long is the answer when you multiply it by a vector?

$4 \times 5 = 20$ entries. The vector must have 5 numbers (one per column), and the answer has 4 (one per row).

Steady state: every tank stays level

Picture each metabolite as a tank of water, and each reaction as a pipe that drains some tanks and fills others. If any tank keeps rising, that molecule piles up until it poisons the cell; if any keeps falling, it runs out and the pathway stops. A healthy, steadily growing cell keeps every tank roughly level. That is a steady state, and in the matrix language it is one short equation:

$$S \cdot v = 0$$

It says: every row of $S$, dotted with the flux vector, comes out zero. Every molecule is made exactly as fast as it is used.

Usually there are many flux vectors that do this. Pyruvate, for example, can go to lactate or to acetate, and any split balances its tank, as long as the two together take everything glycolysis makes. The network has more reactions than balance conditions, so the conditions leave choices open.

Some choices are ruled out by bounds: the lowest and highest speed each reaction is allowed. A one-way reaction can’t run backwards, so its lowest speed is 0. Glucose can come in at most 7.4 mmol per gram dry weight per hour, and acetate can leave at most 6.9. Both limits are borrowed from Mycoplasma pneumoniae, a relative, because nobody has measured them in syn3A.

Flux Plumber turns this into plumbing. Each tank is a molecule and each pipe a reaction; set the valves until every tank is level. The early levels are small made-up networks; the later ones are syn3A’s real core.

Game · Flux Plumber

Keep every tank level

Stars 0/18

Sound on

1. Keep it level2. Split3. Two for one4. The best plan5. syn3A’s glycolysis6. Knockout

Goal: Make the tank stay level: out = in.

1. Keep it level

A tank is a metabolite: a pool of one kind of molecule in the cell. A pipe is a reaction that adds to the pool or takes from it. A valve sets how fast a pipe flows.

Molecules pour in at 5 per second (locked). Set the outflow so the level stops moving. When everything that comes in also goes out, the tank is at steady state.

Start

🔒 In5Out2

Tap a valve (the little wheels) or a chip above to set it.

Check

Ledger: one row per tank, made − used

1. tank X: +1·in − 1·out = 5 − 2 = +3

A made-up network for practice. Flows are in molecules per second; the tank levels show what would happen if the plan ran for a while.

Flux balance analysis

Steady state leaves many plans open, and bounds narrow them down. To pick one, add a goal: make as much of something as possible. For a growing cell, the usual goal is growth itself. In our core, we ask for as much ATP as possible, through the ATPM reaction. Then hand the whole thing to an optimizer. That is flux balance analysis (FBA):

1. Steady state: $S \cdot v = 0$.
2. Bounds: every reaction’s speed between its lowest and highest allowed value.
3. Goal: make one reaction, or a mix of them, as fast as possible.

FBA needs no enzyme speeds and no concentrations, only the reactions and their limits. That is why it can be run on a whole cell’s metabolism.

On syn3A’s core, FBA’s answer is an ATP flux of 21.7 mmol per gram dry weight per hour, about $21.7 \times 1{,}706 \approx 37{,}000$ ATP per second in one cell. On the full model with growth as the goal, it predicts that syn3A doubles every 2.02 hours. Measured syn3A doubles in about 2 hours.

Before the math, look at the picture. Every point on this chart is a plan: how much glucose to take in, and how much acetate to send out. The shaded region holds every plan the limits allow. Press Sweep the goal line and watch the line of equal ATP slide until it touches the last allowed point.

Finding the best plan

Sound on

syn3A’s core with only two choices left: how fast to take up glucose (g) and how fast to send out acetate (a), both in mmol per gram dry weight per hour. Every other flux follows: lactate = 2g − a, oxygen used = a, ATP = 2g + a.

Goal

ATPLactateSave sugar

Knock out

noneackAldh

Sweep the goal line ›

Goal: make the most ATP: 2g + a. Best: the corner (7.4, 6.9), ATP 21.7.

Tap or drag on the chart to place your plan. Hatched means not allowed.

Your plan

glucose g

5

acetate a

3

lactate 2g − a

7

O₂ used = a

3

ATP 2g + a

13

Allowed. About 22,200 ATP per second in one cell.

Glucose limit7.4Acetate limit6.9Oxygen limit15 (no limit)syn3A’s limits (7.4 and 6.9)

Best plan: the corner (7.4, 6.9). All the glucose it can get, as much acetate as can leave, and the rest (7.9) as lactate. ATP = 2 × 7.4 + 6.9 = 21.7 mmol per gram dry weight per hour: about 37,000 ATP per second in one cell.

Every corner of the allowed region

corner (g, a)	lactate	ATP
(7.4, 6.9)	7.9	21.7
(7.4, 0)	14.8	14.8
(3.45, 6.9)	0	13.8
(0, 0)	0	0

In the real model there are 338 directions instead of 2, but the answer still sits on a corner.

Derived from syn3A’s 19-reaction core of the metabolic model iMB155 (Breuer et al. 2019, eLife 8:e36842, CC BY 4.0). The glucose (7.4) and acetate (6.9) limits are borrowed from Mycoplasma pneumoniae (Wodke et al. 2013). The oxygen limit is for exploring: the model sets none. × 1,706 turns mmol per gram dry weight per hour into molecules per second in one syn3A cell.

Math toolkit: linear programming

This toolkit also starts from zero. A linear program is a problem of the form “choose some numbers to make a goal as big as possible, without breaking any limits”, where the goal and the limits are all linear: each number counts with a fixed weight, added up. No squares, no products of two choices. So $2g + a$ is linear, but $g \times a$ or $g^2$ is not. Doubling a choice doubles its effect.

syn3A’s core boils down to two choices. (The other reactions follow from them, through the steady-state rows.) Call them

• $g$: how fast glucose comes in;
• $a$: how fast acetate goes out.

Everything is in mmol per gram dry weight per hour.

The goal. Each glucose sent to lactate makes 2 ATP. Each acetate makes 1 more on top. So

$$\text{ATP} = 2g + a$$

For example, $g = 5$ and $a = 2$ give $2 \times 5 + 2 = 12$.

The limits. Each glucose makes 2 pyruvates. Every pyruvate that doesn’t leave as acetate becomes lactate, so lactate is $2g - a$, and lactate can’t be negative. Together with the two borrowed limits:

1. $g \le 7.4$: glucose can’t come in faster than 7.4.
2. $a \le 6.9$: acetate can’t leave faster than 6.9.
3. $a \le 2g$: you can’t make more acetate than you have pyruvate.
4. $g \ge 0$ and $a \ge 0$: no running backwards.

Now draw it. Put $g$ along the bottom and $a$ up the side. Each limit is a straight line, with the allowed side on one side of it:

• $g = 7.4$ is a vertical line; allowed plans are to its left.
• $a = 6.9$ is a horizontal line; allowed plans are below it.
• $a = 2g$ is a sloping line through $(0, 0)$, $(1, 2)$, $(2, 4)$ and so on; allowed plans are below it, on the side with lots of glucose.

The plans that obey every limit form a region with straight sides, a polygon, called the feasible region. Its corners are where two limits meet. Here there are four: $(0, 0)$; $(7.4, 0)$; $(7.4, 6.9)$; and the point where $a = 6.9$ meets $a = 2g$, which is $g = 6.9 \div 2 = 3.45$, so $(3.45, 6.9)$.

Level lines. Now draw the goal. All the plans that give the same ATP lie on one straight line. For ATP = 10, $2g + a = 10$ passes through $(5, 0)$ and $(2, 6)$, because $2 \times 5 + 0 = 10$ and $2 \times 2 + 6 = 10$. For ATP = 14, it passes through $(7, 0)$ and $(4, 6)$. Bigger ATP means the same slope, shifted up and to the right. Finding the best plan means sliding that line as far up and right as it can go while it still touches the feasible region.

Check yourselfIs the plan g = 3, a = 7 allowed?

No, for two reasons: $a = 7$ is above the acetate limit of 6.9, and it needs more pyruvate than $2 \times 3 = 6$.

Check yourselfHow much ATP does the plan g = 6, a = 4 make, and how much lactate?

ATP: $2 \times 6 + 4 = 16$. Lactate: $2 \times 6 - 4 = 8$.

Why the answer sits on a corner

Here is the feasible region’s four corners, with the ATP each gives:

Corner $(g, a)$	ATP $= 2g + a$
$(0, 0)$	0
$(7.4, 0)$	14.8
$(3.45, 6.9)$	13.8
$(7.4, 6.9)$	21.7

The best is $(7.4, 6.9)$, with 21.7: exactly what FBA found on the full 19-reaction core.

Why must the best plan be a corner? Take any plan inside the region. Moving it a little in the direction of more glucose raises ATP, and nothing stops you until you hit a limit, an edge. Now walk along that edge. The goal is a straight line, so as you walk, ATP changes at a steady rate: always up, always down, or not at all. If it goes up, keep walking until the next limit, which is a corner. So the best plan always sits on a corner, where two limits meet. A solver doesn’t need to try every plan, only the corners, and clever ones skip most of those too.

Read the winning corner in words: take in all the glucose the cell can get (7.4), send out as much acetate as can leave (6.9), and turn the rest of the pyruvate into lactate ($2 \times 7.4 - 6.9 = 7.9$). ATP: $2 \times 7.4 + 6.9 = 21.7$, or $21.7 \div 7.4 \approx 2.9$ per glucose. The limits that the answer presses against, here the glucose and acetate limits, are called binding constraints. They tell you what would help: raise the acetate limit to 10 and ATP rises to $14.8 + 10 = 24.8$.

Sometimes there is a tie. If the goal’s level lines are parallel to one edge of the region, every point along that edge is equally good: an alternative optimum. The full syn3A model has one. It lets pyruvate itself leak out of the cell, and if pyruvate leaks instead of becoming lactate, the extra NADH can be emptied by the oxidase. That gives exactly the same ATP. The solver has to pick one, and it happens to pick the leak, so the full model’s FBA shows no lactate at all. Whenever an FBA answer looks odd, check for a tie.

Now run FBA yourself. The starter builds the matrix; you set up the goal and the bounds, call the solver, and return the answer. The checks expect 21.7 for syn3A’s core, and 14.8 with the acetate exit closed.

Python exercise

Flux balance analysis on syn3A

from virtualcell.constants import FLUX_TO_PER_CELL_PER_SECOND # 1,706 molecules per second from virtualcell.metabolic_model import core_reactions, solve_lp def flux_balance(reactions, objective, changes=None): """Flux balance analysis: the steady-state fluxes that make `objective` as fast as possible.""" # Step 1: metabolites, ids, S = stoichiometric_matrix(reactions) # Step 2: the objective as a list c, one number per reaction: 1 for the # reaction whose id is objective, 0 for every other one. # Step 3: the bounds, a list of (lb, ub) pairs, one per reaction. If changes # has an entry for a reaction's id, use that pair instead, for example # changes = {"ACout": (0, 0)} closes the acetate exit. changes may be # None: then nothing changes. # Step 4: steady state means S·v = 0, so b_eq is a list of zeros, one per # metabolite. Then: result = solve_lp(c, S, b_eq, bounds) # (solve_lp makes c·v as large as possible: it maximizes.) # Step 5: return {"status": result["status"], "objective": result["objective"], # "fluxes": a dict from each reaction id to its flux}. # The fluxes are result["x"], in the same order as ids. When the status # isn't "optimal", result["x"] is None: then return None for "fluxes". raise NotImplementedError("Write your code here") def stoichiometric_matrix(reactions): """The same as your answer to exercise 5.""" metabolites = sorted({name for reaction in reactions for name in reaction["stoich"]}) reaction_ids = [reaction["id"] for reaction in reactions] S = [[reaction["stoich"].get(name, 0) for reaction in reactions] for name in metabolites] return metabolites, reaction_ids, S # A small made-up network: food A comes in (at most 10), becomes B, and B splits # into C (what the cell wants) and a byproduct D, which can leave at most 4. TOY = [ {"id": "R1", "stoich": {"A": 1}, "lb": 0, "ub": 10}, {"id": "R2", "stoich": {"A": -1, "B": 1}, "lb": 0, "ub": 1000}, {"id": "R3", "stoich": {"B": -1, "C": 1, "D": 1}, "lb": 0, "ub": 1000}, {"id": "R4", "stoich": {"C": -1}, "lb": 0, "ub": 1000}, {"id": "R5", "stoich": {"D": -1}, "lb": 0, "ub": 4}, ] # Remove the #s below to try it: first the toy, then syn3A's core. # print(flux_balance(TOY, "R4")) # result = flux_balance(core_reactions(), "ATPM") # print("Fluxes in mmol per gram dry weight per hour (molecules per second in one cell):") # for rid, flux in result["fluxes"].items(): # print(f"{rid:>7}: {flux:6.2f} ({flux * FLUX_TO_PER_CELL_PER_SECOND:6,.0f} per s)") # print("Without the acetate exit:", flux_balance(core_reactions(), "ATPM", {"ACout": (0, 0)})["objective"])

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfIf acetate couldn't leave at all (acetate limit 0), what would the best plan be?

The region shrinks to the line $a = 0$ from $g = 0$ to $g = 7.4$. The best corner is $(7.4, 0)$: all lactate, $2 \times 7.4 = 14.8$ ATP.

Check yourselfWhy is the corner (3.45, 6.9) where it is?

It is where the acetate limit, $a = 6.9$, meets the pyruvate limit, $a = 2g$. Solving $6.9 = 2g$ gives $g = 3.45$.

Knockouts: which genes can syn3A lose

FBA can also ask: what if a gene were missing? Each reaction in the model has a gene rule saying which genes it needs. “A and B” means it needs both, like the PTS, which needs all four of ptsI, ptsH, crr and ptsG. “A or B” means either will do, for two enzymes that can stand in for each other. A reaction with no gene (such as the oxidase) is assumed to always work.

To knock out a gene in the computer, find every reaction whose rule fails without it, set that reaction’s speed to zero, and solve again. If the best goal drops to zero (or close to it), the gene is essential: the cell can’t live without it. This kind of test, done in the computer rather than in a dish, is called in silico.

On syn3A’s core:

• Knock out any of the 13 genes of the PTS and glycolysis, and ATP drops to 0. There is no other way for sugar to get through.
• Knock out pdhC, lpdA, pta or ackA, and the acetate route closes. ATP drops to 14.8 (32% less): all lactate.
• Knock out ldh, and pyruvate has nowhere to go but acetate. Acetate can leave at only 6.9, so glucose has to drop to 3.45, and ATP falls to 13.8 (36% less).

On the full model with growth as the goal, 123 of its 155 genes are essential: 79%. The acetate genes aren’t essential, but losing one slows growth from a doubling time of 2.02 hours to 3.22 hours.

How do you check that in real cells? The team at the J. Craig Venter Institute used transposon mutagenesis. A transposon is a piece of DNA that can insert itself at random places in a genome; wherever it lands inside a gene, it breaks that gene. Make a large population of cells with insertions in random places, let them grow, and read where the insertions are. Genes that never carry an insertion in growing cells are essential. Genes whose broken cells grow, but poorly, are quasi-essential. The rest are non-essential.

For 153 of the model’s genes, here is how the model and the lab compare, counting quasi-essential genes as essential:

	Lab: essential or quasi	Lab: non-essential
Model: essential	123	0
Model: non-essential	18	12

The model is right for $123 + 12 = 135$ genes: an accuracy of $135 \div 153 = 88\%$. Of the 141 genes the lab says the cell needs, it finds 123 (87%, its sensitivity). Of the 12 the lab says it can lose, it gets all 12 (100%, its specificity). Every gene the model calls essential is at least quasi-essential in the lab.

Accuracy can flatter a model when one answer is much more common than the other. A model that called every gene essential would score $141 \div 153 = 92\%$ here and know nothing. So modelers also use the Matthews correlation coefficient (MCC), a score from −1 to 1 that is 1 for a perfect prediction, 0 for no better than guessing, and −1 for always wrong. syn3A’s model scores about 0.59: clearly better than guessing, clearly not perfect.

The misses are the interesting part. The lab finds four genes essential that the model thinks the cell can lose, and all four sit right next to glycolysis:

• ldh (0475), lactate dehydrogenase. Without it, the model lets pyruvate leak out of the cell and empties the extra NADH with the oxidase, the one with no known gene. Real cells can’t do that, it seems.
• pta (0229) and ackA (0230), the acetate route. The model says the cell would only grow slower, 3.22 hours instead of 2.02. The real cells die.
• gapdh (0451), a second GAPDH that makes NADPH instead of NADH. Breuer and colleagues suggest the cell needs NADPH in ways the model doesn’t capture.

Each miss points to something we don’t yet understand about syn3A.

This last exercise is a challenge. Write the gene-rule check, the knockout, and the scan over every gene, then find syn3A’s essential core genes. The success message compares your list with the lab’s.

Python exercise

Challenge: which genes can syn3A lose?

from virtualcell.metabolic_model import GENES, core_reactions, solve_lp def reaction_allowed(rule, knocked_out): """Can a reaction still run when the genes in knocked_out are gone?""" # A rule is made of gene names, "and", "or" and brackets, such as "(A and B) or C": # "and" means both are needed (two parts of one enzyme), "or" means either # will do (two enzymes that do the same job). # Python can work out a rule like that once every gene name is replaced by # True (still there) or False (knocked out): (True and False) or True is True. # Step 1: an empty rule (no known gene) is always allowed: return True. # Step 2: put spaces around the brackets, so that split() keeps them apart: # words = rule.replace("(", " ( ").replace(")", " ) ").split() # Step 3: make a new list of words: "and", "or", "(" and ")" stay as they are; # a gene becomes "False" if it is in knocked_out, and "True" if not. # Step 4: join the words with spaces and let Python work it out: # return eval(" ".join(new_words)) # (eval runs a string as Python. It is safe here, because the string # holds nothing but True, False, and, or and brackets.) raise NotImplementedError("Write your code here") def knockout(reactions, genes): """The reactions with these genes removed, as a new list (the input is not changed).""" # For each reaction, make a copy: dict(reaction). If reaction_allowed(its # "rule", genes) is False, set the copy's "lb" and "ub" to 0. # Return the list of copies. raise NotImplementedError("Write your code here") def essential_genes(reactions, objective, fraction=0.1): """The genes the cell can't live without: a sorted list of locus tags.""" # Step 1: the normal best: flux_balance(reactions, objective)["objective"]. # Step 2: every gene, once: the "genes" lists of all the reactions, in a set. # Step 3: for each gene, knock it out alone and run flux_balance again. The # gene is essential if the status isn't "optimal", or if the best # objective falls below fraction × the normal best. # Return the essential genes, sorted. raise NotImplementedError("Write your code here") def flux_balance(reactions, objective, changes=None): """The same as your answer to exercise 7.""" changes = changes or {} metabolites, ids, S = stoichiometric_matrix(reactions) c = [1 if rid == objective else 0 for rid in ids] bounds = [changes.get(r["id"], (r["lb"], r["ub"])) for r in reactions] result = solve_lp(c, S, [0] * len(metabolites), bounds) fluxes = dict(zip(ids, result["x"])) if result["status"] == "optimal" else None return {"status": result["status"], "objective": result["objective"], "fluxes": fluxes} def stoichiometric_matrix(reactions): """The same as your answer to exercise 5.""" metabolites = sorted({name for reaction in reactions for name in reaction["stoich"]}) reaction_ids = [reaction["id"] for reaction in reactions] S = [[reaction["stoich"].get(name, 0) for reaction in reactions] for name in metabolites] return metabolites, reaction_ids, S # Remove the #s below to try it on syn3A's core. GENES has each gene's name. # print(reaction_allowed("JCVISYN3A_0227 and JCVISYN3A_0228", {"JCVISYN3A_0227"})) # False # core = core_reactions() # for tag in essential_genes(core, "ATPM"): # print(tag, GENES[tag]["name"], "is essential") # Optional extra: the whole published model, iMB155 (338 reactions, 155 genes). # It is far too big for solve_lp in plain Python, so this uses scipy's solver. # Warning: the first time, the browser downloads scipy, about 17 MB; on a phone # or a slow connection, run it on your computer instead. Remove the #s to try it. # import scipy # from virtualcell.metabolic_model import IMB155_OBJECTIVE, imb155_reactions, solve_lp_scipy # solve_lp = solve_lp_scipy # from now on, flux_balance above uses scipy # full = imb155_reactions() # growth = flux_balance(full, IMB155_OBJECTIVE)["objective"] # print(f"growth rate {growth:.3f} per hour: one doubling every {0.6931 / growth:.2f} hours") # lost = essential_genes(full, IMB155_OBJECTIVE) # 156 linear programs: about 5 seconds # print(len(lost), "of 155 genes are essential by plain knockout") # # Some reactions must carry a minimum flux (lb > 0 or ub < 0), such as the ATP # # synthase pumping protons. Switching one off breaks that rule, so count it as lethal: # must = {g for r in full if r["lb"] > 0 or r["ub"] < 0 for g in r["genes"]} # print(len(set(lost) | must), "of 155 counting those: the paper says 123 (79%), and 88% agree with the lab")

Ctrl+Enter runs your code. Tab indents; press Esc then Tab to leave the editor.

Run and checkHint (3 left)Show solutionReset

Check yourselfA reaction's rule is “A and B”. Can it run if B is knocked out? What if the rule is “A or B”?

With “A and B”, no: it needs both. With “A or B”, yes: A can do the job alone.

Check yourselfA model calls 50 genes essential, and the lab agrees for 45 of them. It calls 50 non-essential, and the lab agrees for 40. What is its accuracy?

$(45 + 40) \div 100 = 85\%$.

Code milestone: flux balance analysis

Two new modules join the engine in this lesson:

• virtualcell/enzymes.py holds michaelis_menten, max_rate, busy_fraction, competitive_rate, noncompetitive_rate and fit_michaelis_menten.
• virtualcell/metabolism.py holds stoichiometric_matrix, mass_balance, is_steady_state, flux_balance, reaction_allowed, knockout and essential_genes.

One is written for you. virtualcell/metabolic_model.py has the real core of syn3A’s metabolism (19 reactions from iMB155, Breuer et al. 2019, CC BY 4.0), the five enzymes’ constants, a linear-program solver in plain Python, and a loader for the whole published model. constants.py has this lesson’s numbers, including the 1,706 molecules per second per unit of flux.

Copy each function you wrote into its module, replacing the stub, then run the checks from the engine folder:

python -m pytest

When they pass, the engine can run syn3A’s metabolism:

from virtualcell.enzymes import michaelis_menten
from virtualcell.metabolic_model import core_reactions
from virtualcell.metabolism import flux_balance, essential_genes

print(michaelis_menten(0.28, 929, 0.28))      # 464.5, half of PGI's top speed

core = core_reactions()
best = flux_balance(core, "ATPM")
print(best["objective"])                      # 21.7
print(best["fluxes"]["LACout"])               # 7.9

closed = flux_balance(core, "ATPM", {"ACout": (0, 0)})
print(closed["objective"])                    # 14.8, no acetate

print(len(essential_genes(core, "ATPM")))     # 13

For the whole model, 338 reactions, metabolic_model.solve_lp_scipy does the same job with SciPy’s faster solver; the knockouts exercise has an optional extra that uses it.

Numbers to remember

Quantity	Value
syn3A’s metabolic model (iMB155)	338 reactions, 304 metabolites, 155 genes
Typical enzyme	$k_\text{cat}$ about 10 per second, $K_\text{m}$ about 0.1 mM
Record speed-up	$10^{17}$ (78 million years → about 25 ms)
Barrier and speed	5.9 kJ/mol lower = 10× faster at 37 °C
Michaelis–Menten	$v = V_\text{max}[\text{S}]/(K_\text{m} + [\text{S}])$; half speed at $[\text{S}] = K_\text{m}$, 90% at $9K_\text{m}$
Competitive / noncompetitive	raises $K_\text{m}$ / lowers $V_\text{max}$
Glycolysis	2 ATP and 2 NADH per glucose; up to 4 ATP with acetate
syn3A’s glucose	about 15,000 molecules per second, through the PTS (no hexokinase)
NAD⁺ without recycling	runs out in about 1.5 seconds
ATP	about 50 kJ/mol, about 20 kT; about 74,000 in a syn3A cell, remade every 2–3 s
1 mM in syn3A	about 20,000 molecules
1 mmol per gram dry weight per hour	about 1,706 molecules per second in one syn3A cell
Steady state	$S \cdot v = 0$
FBA on syn3A’s core	21.7 ATP per gram dry weight per hour (about 37,000 per second)
Gene essentiality	123 of 155 essential in the model; 88% agreement with the lab

Next lesson: gene regulation and circuits. You will see how genes switch each other on and off, build switches and clocks out of feedback, and write an engine that simulates them.

Sources and further reading

• Breuer M, et al. Essential metabolism for a minimal cell. eLife 8:e36842 (2019). The metabolic model iMB155, its FBA results and gene essentiality; CC BY 4.0. Supplementary files 3 (essentiality), 7 (kcat and protein copies) and 10 (the model).
• Thornburg ZR, et al. Fundamental behaviors emerge from simulations of a living minimal cell. Cell 185:345–360 (2022). Glucose uptake, ATP production and costs, starting concentrations, and the $K_\text{m}$ values.
• Thornburg ZR, et al. Kinetic modeling of the genetic information processes in a minimal cell. Frontiers in Molecular Biosciences 6:130 (2019). The ATP cost of gene expression.
• Hutchison CA III, et al. Design and synthesis of a minimal bacterial genome. Science 351:aad6253 (2016). The transposon experiments behind the essential-gene calls.
• Wodke JAH, et al. Dissecting the energy metabolism in Mycoplasma pneumoniae through genome-scale metabolic modeling. Molecular Systems Biology 9:653 (2013). The borrowed glucose and acetate limits.
• Bar-Even A, et al. The moderately efficient enzyme: evolutionary and physicochemical trends shaping enzyme parameters. Biochemistry 50:4402–4410 (2011).
• Radzicka A, Wolfenden R. A proficient enzyme. Science 267:90–93 (1995). The 78-million-year reaction.
• Michaelis L, Menten ML. Die Kinetik der Invertinwirkung. Biochemische Zeitschrift 49:333–369 (1913).
• Johnson KA, Goody RS. The original Michaelis constant: translation of the 1913 Michaelis–Menten paper. Biochemistry 50:8264–8269 (2011).
• Bennett BD, et al. Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nature Chemical Biology 5:593–599 (2009).
• Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nature Biotechnology 28:245–248 (2010). A short, friendly introduction.
• syn3A’s genome and genes: NCBI GenBank record CP016816.2.
• Milo R, Phillips R. Cell Biology by the Numbers. Garland Science, 2015. Free online at book.bionumbers.org. Enzyme rates, reactions without enzymes, and the energy in ATP.