Skip to content

A physiologically-based digital twin for alcohol consumption—predicting real-life drinking responses and long-term plasma PEth

This is an in house modelling project where a model governing the dynamics of alcohol rate of appearance in plasma, following consumption of different alcoholic beverages paired with either:

  • i) non-alcoholic beverages
  • ii) food
  • iii) other alcoholic beverages
  • iV) solely the alcoholic drink itself

The article is published in NPJ digital twin, Podéus, H., Simonsson, C., Nasr, P. et al. A physiologically-based digital twin for alcohol consumption—predicting real-life drinking responses and long-term plasma PEth. npj Digit. Med. 7, 112 (2024), and can be found following this DOI: https://doi.org/10.1038/s41746-024-01089-6.

The model

model figure

This model describe how beverages pass the stomach, and if they have an alcoholic content, how the alcohol is taken up into the intestines, transported into the blood, and then eliminated via enzymes in the liver. The rate of gastric emptying is influenced by the amount of available caloric content in the stomach and if there is food present.

The model is govern by multiple type of inputs.

  • There are anthropometric constants (sex, height, and weight) set by a "constant" in the activity object.
  • The food and drink variables (alcohol concentration, rate of consumption, caloric content, and consumption duration) is govern by "piecewise_constant" in the activity object, as it is intent for multiple drinks and meals to be possible to occur in sequence.

To correctly keep track of the digestion of food and beverage kcal in the stomach, the model also includes events to: update the state goering the rate of gastric emptying with the remaining kcal (that has not been digested) once a new drink starts, and an event that similar updates the solid kcal, from food, once a new meal is consumed - as the digestion is dependent on the starting amount of kcal.

The model file
########## NAME
alcoholmodel
########## METADATA
time_unit = m
########## MACROS
########## STATES
d/dt(Vol_Stomach) = + vol_drink_per_time*10 - r2            // dL
d/dt(Kcal_Liquid) = + vol_drink_per_time*kcal_liquid_per_vol - rKcal_clearance         // kcal [ L * kcal/L = kcal]

d/dt(max_Kcal_Solid) = 0
d/dt(Kcal_Solid)     = r_Kcal_Solid

d/dt(EtOH_Pool) = + r_poolIn*Vol_Stomach - r_poolOut                    //mg

d/dt(Conc_EtOH_Stomach) = r_drinkEtOH  - r_poolIn + r_poolOut/Vol_Stomach // mg/dL
d/dt(Mass_EtOH_Intestines) = + r2*Conc_EtOH_Stomach - r3  - r4          // mg

d/dt(Blood_Conc) =  + (r3/V_Blood) - r5*(V_Liver/V_Blood)               //mg/dL 

d/dt(Plasma_acetate) = + r5 - r6                  //mg/dL  

d/dt(PEth)       = rPEth - rPEth_bound + rPEth_release - rPEth_clearance 
d/dt(PEth_Bound) = rPEth_bound - rPEth_release

d/dt(Kcal_remain) = + vol_drink_per_time*kcal_liquid_per_vol - r2*kcal_liquid_per_vol/10          // kcal [ L * kcal/L = kcal]

d/dt(time_elapsed) = 1

Vol_Stomach(0)          = 0.001
Kcal_Liquid(0)          = 0
max_Kcal_Solid(0)       = 0
Kcal_Solid(0)           = 0
EtOH_Pool(0)            = 0
Conc_EtOH_Stomach(0)    = 0
Mass_EtOH_Intestines(0) = 0
Blood_Conc(0)           = 0
Plasma_acetate(0)       = 0
PEth(0)                 = 0
PEth_Bound(0)           = 0
Kcal_remain(0)          = 0
time_elapsed(0)         = 0

########## PARAMETERS
kPEth            = 14824.3201129844
kPEth_out        = 155.700230807449
kPEth_bind       = 16899.8070912636
kPEth_release    = 0.00567465407687899
k_poolIn         = 8.51709430245296
k_poolOut        = 0.346199847345036
Vmax             = 2018.19466075863
km               = 16336.3631792709
k_kcal           = 154.363454312840
k3               = 84.3794789018948
k4               = 852.816905367382
k6               = 0.127255213814248
VmaxADH          = 0.935326859734748
VmaxCYP2E1       = 0.177371543132190
KmADH            = 8.97358830896085  
KmCYP2E1         = 40.9840079864542  
k_Kcal_clearance = 0.00570193833670537

########## VARIABLES
SS_vol = 0.001

V_Liver = 15
V_Blood = ((1-sex)*(0.3561 * height^3 + 0.03308 * weight + 0.1833) + sex*(0.3669 * height^3 + 0.03219 * weight + 0.6041))*10 // dL     Nadler's Equation for total blood volume

conc_drink = (EtOH_conc*789.1) //mg/dL

vADH    = (VmaxADH*Blood_Conc/( KmADH + Blood_Conc))           
vCYP2E1 = (VmaxCYP2E1*Blood_Conc/( KmCYP2E1 + Blood_Conc))  

kcal_solid_vol = max(1,Kcal_Solid)/(4*100) // 4kcal/g assumes 1g/ml, 100 scales ml--> dL 

// Model reactions
r2 = Vmax*( (Vol_Stomach-SS_vol)/((Vol_Stomach-SS_vol) + km) ) * exp( - max(0, Kcal_Liquid)/k_kcal )

r_Kcal_Solid = (max_Kcal_Solid * -1.88 * 0.010* (0.010*max(0, time_elapsed))^0.86  * exp(-(0.010*max(0, time_elapsed))^1.86))    // doi: 10.1016/S0002-9270(00)00868-6

r_poolIn  = gt(Kcal_Solid, 1e-3) * max(0, (Conc_EtOH_Stomach - EtOH_Pool/kcal_solid_vol) * k_poolIn)
r_poolOut =  gt(Kcal_Solid, 1e-3) * EtOH_Pool * k_poolOut

r_drinkEtOH     = (conc_drink - Conc_EtOH_Stomach)*(vol_drink_per_time*10/Vol_Stomach)
rKcal_clearance = Kcal_Liquid*k_Kcal_clearance

r3 = Mass_EtOH_Intestines*k3 
r4 = Mass_EtOH_Intestines*k4

r5  = vADH + vCYP2E1

r6 = Plasma_acetate*k6

rPEth           = kPEth*Blood_Conc
rPEth_clearance = kPEth_out*PEth
rPEth_bound     = kPEth_bind*PEth
rPEth_release = max(0, kPEth_release*( PEth_Bound - PEth ))

########## FUNCTIONS

########## EVENTS
new_drink = gt(vol_drink_per_time,0), Kcal_Liquid, Kcal_remain, Kcal_remain, 0
new_meal  = gt(kcal_solid,0), max_Kcal_Solid, kcal_solid+Kcal_Solid, Kcal_Solid, kcal_solid+Kcal_Solid, time_elapsed, 0 

########## OUTPUTS

########## INPUTS
EtOH_conc           = EtOH_conc           @ 0
vol_drink_per_time  = vol_drink_per_time  @ 0
kcal_liquid_per_vol = kcal_liquid_per_vol @ 0
kcal_solid          = kcal_solid          @ 0
drink_length        = drink_length        @ 0

sex                 = sex                 @ 1       // male=1 , female=0
weight              = weight              @ 104     // Kg
height              = height              @ 1.85    // m

########## FEATURES

Blood alcohol concentration = Blood_Conc    [mg/dL]
Acetate in plasma = Plasma_acetate/10.2     [mM]
Breath alcohol concentration (g/210L) = 0.840*(Blood_Conc/1000) + 0.00367 [g/210L]
Blood BrAC concentration (g/dL) = Blood_Conc/1000 [g/dL]

Gastric volume = Vol_Stomach*100 [mL]
PEth = PEth [ng/mL]

Implementation

An python implementation of the model, including a simulation with a) an alcoholic drink, b) an alcoholic drink paired with a meal, c) an non-alcoholic drink, and d) an alcoholic drink and plotting the corresponding results for the observables a) blood alcohol concentration (BAC), b) blood alcohol concentration, c) gastric volume, and d) phosphatidylethanol (PEth) is given below.

The python implementation is shown below and can be downloaded here.

Requirements

The python code requires the following packages to be installed.

uv add numpy matplotlib
pip install numpy matplotlib

Example code

import matplotlib.pyplot as plt
import numpy

import sund


#%% Install the model by using the SUND function 'install_model'
MODEL_NAME = "alcoholmodel"
sund.install_model(f"{MODEL_NAME}.txt")

#%% Load the model by using the sund function 'load_model'
model = sund.load_model(f"{MODEL_NAME}")

#%% Create an activity and simulate an alcoholic drink
activity = sund.Activity(time_unit=model.time_unit)
activity.add_output("piecewise_constant", "EtOH_conc",           t = [0, 20], f=[0, 5.1, 0])
activity.add_output("piecewise_constant", "vol_drink_per_time",  t = [0, 20], f=[0, 0.05135, 0])
activity.add_output("piecewise_constant", "kcal_liquid_per_vol", t = [0, 20], f=[0, 129.82, 0])
activity.add_output("piecewise_constant", "drink_length",        t = [0, 20], f=[0, 20, 0])
activity.add_output("constant", "sex",    f=1)
activity.add_output("constant", "weight", f=82.66)
activity.add_output("constant", "height", f=1.77)

# Create a simulation object
sim = sund.Simulation(model, activities = activity, time_unit = model.time_unit, time_vector = numpy.arange(0, 250, 0.1))

sim.simulate()

# Plot the results
fig = plt.figure(figsize=(18*(1/2.54), 18*(1/2.54)))
ax = fig.add_subplot(2, 2, 1)
ax.plot(sim.time_vector, sim.feature_data[:, sim.feature_names.index('Blood alcohol concentration')], label='Blood alcohol concentration')
ax.set_xlabel('Time (min)')
ax.set_ylabel('BAC (mg/dL)')
ax.legend(edgecolor='None', loc='upper right')

#%% Create an activity and simulate an alcoholic drink with food
activity = sund.Activity(time_unit=model.time_unit)
activity.add_output("piecewise_constant", "EtOH_conc",           t = [20, 35],   f=[0, 20, 0])
activity.add_output("piecewise_constant", "vol_drink_per_time",  t = [20, 35],   f=[0, 0.009606, 0 ])
activity.add_output("piecewise_constant", "kcal_liquid_per_vol", t = [20, 35],   f=[0, 360, 0])
activity.add_output("piecewise_constant", "drink_length",        t = [20, 35],   f=[0, 15, 0])
activity.add_output("piecewise_constant", "kcal_solid",          t = [1e-8, 20], f=[0, 700, 0])
activity.add_output("constant", "sex",    f=1)
activity.add_output("constant", "weight", f=75.80)
activity.add_output("constant", "height", f=1.83)

# Create a simulation object
sim = sund.Simulation(model, activities = activity, time_unit = model.time_unit, time_vector = numpy.arange(0, 400, 0.1))

sim.simulate()

# Plot the results
ax = fig.add_subplot(2, 2, 2)
ax.plot(sim.time_vector, sim.feature_data[:, sim.feature_names.index('Blood alcohol concentration')], label='Blood alcohol concentration')
ax.set_xlabel('Time (min)')
ax.set_ylabel('BAC (mg/dL)')
ax.legend(edgecolor='None', loc='upper right')

#%% Create an activity and simulate a non-alcoholic drink
activity = sund.Activity(time_unit=model.time_unit)
activity.add_output("piecewise_constant", "vol_drink_per_time",  t = [0, 3], f=[0, 0.5/3, 0])
activity.add_output("piecewise_constant", "kcal_liquid_per_vol", t = [0, 3], f=[0, 664, 0])
activity.add_output("piecewise_constant", "drink_length",        t = [0, 3], f=[0, 3, 0])
activity.add_output("constant", "sex",    f=1)
activity.add_output("constant", "weight", f=65.00)
activity.add_output("constant", "height", f=1.71)

# Create a simulation object
sim = sund.Simulation(model, activities = activity, time_unit = model.time_unit, time_vector = numpy.arange(0, 125, 0.1))

sim.simulate()

# Plot the results
ax = fig.add_subplot(2, 2, 3)
ax.plot(sim.time_vector, sim.feature_data[:, sim.feature_names.index('Gastric volume')], label='Gastric volume')
ax.set_xlabel('Time (min)')
ax.set_ylabel('Volume (mL)')
ax.legend(edgecolor='None', loc='upper right')

#%% Create an activity and simulate the PEth profile from an alcoholic drink
activity = sund.Activity(time_unit=model.time_unit)
activity.add_output("piecewise_constant", "EtOH_conc",           t = [0, 15], f=[0, 6.41, 0])
activity.add_output("piecewise_constant", "vol_drink_per_time",  t = [0, 15], f=[0, 0.0473173, 0])
activity.add_output("piecewise_constant", "kcal_liquid_per_vol", t = [0, 15], f=[0, 322.3, 0])
activity.add_output("piecewise_constant", "drink_length",        t = [0, 15], f=[0, 15, 0])
activity.add_output("constant", "sex",    f=1)
activity.add_output("constant", "weight", f=59.23)
activity.add_output("constant", "height", f=1.78)

# Create a simulation object
sim = sund.Simulation(model, activities = activity, time_unit = model.time_unit, time_vector = numpy.arange(0, 400, 0.1))

sim.simulate()

# Plot the results
ax = fig.add_subplot(2, 2, 4)
ax.plot(sim.time_vector, sim.feature_data[:, sim.feature_names.index('PEth')], label='PEth')
ax.set_xlabel('Time (min)')
ax.set_ylabel('PEth (ng/mL)')
ax.legend(edgecolor='None', loc='upper right')


plt.show()