Create Aluminum MRLs in Matlab

One of the most cited papers defending the safety and use of Aluminum (Al) in vaccines was published by [Keith in 2002]. They used a mathematical models and various assumptions to model the infant body burdens during the first year of life for breast milk and formula diets and for a standard vaccination schedule. They then compared those body burdens with that expected for intake at a level considered safe for intermediate-duration exposure. The paper properly adjusts dietary exposures and assumes that only 0.78% of dietary aluminum makes it into the bloodstream (the rest is excreted).

In this post, I share with you how to re-create the Minimum Risk Levels (MRL) computed by Keith using Matlab, which I use to validate my understanding of their models and how they arrived at their result. An MRL is the dose which is expected to be safe for human exposure. I was interested in this model because I wanted to better understand the sensitivity of the model to various alternative MRL’s that have been derived. For example, the paper uses 2mg/kg per day of Al as the most that a person can consume without showing any signs of toxicity. However, the European Food Safety Authority has issued a statement saying that at most 1 mg/kg per week (or 0.14 mg/kg per day) is tolerable (Tolerable Weekly Intake, TWI). The latter would suggest that the MRL curves in the paper could potentially be 10 times smaller, but the paper does not offer any kind of sensitivity analysis on this assumption, so being able to model this ourselves provides an additional window of insight. It also isn’t not clear what the Al burden is if you add dietary and vaccine sources together. There are a number of other weaknesses with the paper, but for the moment I’m only interested in the model as if the assumptions are true.

Below is a graph that uses the original graph and then overlays a Matlab-generated graph over the same coordinates. If my understanding of the paper is correct, then the colorful lines which were generated in Matlab should closely match the black-and-white lines that are part of the original paper. The lines fit nearly perfectly, even if my code was based on a discrete, daily accumulation of exposures (vs. continuous as modeled in the paper).

validate_keith

The only discrepancy is that observed for Formula in the second half of the year, which to me indicates that there are some details that weren’t clearly disclosed, but pertained to Formula intake only. I believe this small difference is unlikely to affect our interpretation of the sensitivities to the underlying assumptions.

Exploring these curves for various alternative MRLs using more conservative assumptions will be the subject of a future post.

Matlab Code

First download the Background Image.

Also save the following function into a file called priest.m

function p = priest(t, Q)
    t = t - t(1) + 1;
    p = repmat(Q, length(t), 1) ...
        .* repmat([1; 0.354*t(2:end).^(-0.32)], 1, size(Q, 2));
end

Then run the following script with the above background image in the working directory. What the code does is to model each day’s exposure separately. Priest’s function tells you the percentage of the exposure that is still in your body on consecutive days, so at a high level, what we’re doing here is just adding up each day’s exposure to get a total cumulative amount.

% Load in Keith's graph as an image
close all hidden;
ha = axes('units','normalized', ...
            'position',[0 0 1 1]);
% Move the background axes to the bottom
uistack(ha,'bottom');
% Load in a background image and display it using the correct colors
I=imread('KeithMRL.png');
hi = imagesc(I);
colormap gray
% Turn the handlevisibility off so that we don't inadvertently plot into the axes again
% Also, make the axes invisible
set(ha,'handlevisibility','off', ...
            'visible','off')
% Now we can use the figure, as required.
% Create a graph on top of the same coordinates as used by Keith.
h = axes('position',[0.1423,0.2342,0.5228,0.7208], 'Color', 'None');

% ----------------------------------------------------------------------
% PARAMETERS
% ----------------------------------------------------------------------
noael = 62; % oral daily dose in mg Al
mrl = noael /(3*10);
t = (1:400)'; % time in days following uptake
uf = 0.0078; % uptake factor
bw = 3.23 + 10*(1-exp(-0.0028*t)); % body weight
breast_milk_concentrations = 0.040; % mg per liter
semisolid_food = 0.7; %mg
formula_concentrations = 0.225; % mg per liter
fluid_intake = [linspace(0.67, 0.9, 6*30+1)'; repmat(0.900, ...
    length(t)-6*30-1, 1)]; % liters per day.

% ------------------------------------------------------------------------
% CALCULATE THE VARIOUS EXPOSURES
% ------------------------------------------------------------------------
nmrl = size(mrl, 2);

% Compute the MRL
MRL = zeros(length(t), nmrl);
for i = 1:length(t)
    Q = uf*mrl.*bw(i);
    MRL = MRL + [zeros(i-1,nmrl); priest(t(i:end), Q)];
end

% Table of Al exposure from Vaccines at different months:
Q = [0, 0.25;
    2, 1.2;
    4, 0.85;
    6, 1.2;
    12, 0.85];
V = zeros(length(t), 1);
for iQ = 1:size(Q, 1)
    tQ = Q(iQ,1)*30+1;
    qQ = Q(iQ, 2);
    V = V + [zeros(tQ-1,1); priest(t(tQ:end), qQ)];
end

% breast milk
milk = zeros(length(t), 1);
for i = 1:length(t)
    Q = uf*fluid_intake(i).*breast_milk_concentrations;
    if (i >= 30*6)
        % add semi-solid food in the second half of the year.
        Q = Q + uf*semisolid_food;
    end
    milk = milk + [zeros(i-1,1); priest(t(i:end), Q)];
end

% Formula
frmla = zeros(length(t), 1);
for i = 1:length(t)
    Q = uf*fluid_intake(i).*formula_concentrations;
    if (i >= 30*6)
        % add semi-solid food in the second half of the year.
        Q = Q + uf*semisolid_food;
    end
    frmla = frmla + [zeros(i-1,1); priest(t(i:end), Q)];
end

% ------------------------------------------------------------------------
% CREATE THE GRAPH
% ------------------------------------------------------------------------
semilogy(t, MRL, 'm--',...
                'LineWidth',1.5,...
                'MarkerEdgeColor','m');
hold on;
semilogy(t, V, 'g-',...
    'LineWidth',1.5,...
    'MarkerEdgeColor','g');
semilogy(t, milk, 'b-',...
    'LineWidth',1.5,...
    'MarkerEdgeColor','b');
semilogy(t, frmla, 'b:',...
    'LineWidth',1.5,...
    'MarkerEdgeColor','b');

hold off;
set(h, 'YLim', [10^(-4), 10]);
% Make the graph transparent.
set(h, 'Color', 'None');
legend({'MRL Golub', 'Vaccines', 'Breast Milk', 'Formula'}, 'location', 'southeast');


print(1,'-dpng','validate_keith');

 

[1] Keith, L.S., Jones, Chou; Aluminum Toxicokinetics regarding infant diet and vaccinations; Vaccine 20 (2002) S13-S17 PMID: 12184359  
[2]  EFSA AFC (Former Panel on food additives, flavourings, processing aids and materials in contact with food), 2008. Safety of aluminium from dietary intake[1] – Scientific Opinion of the Panel on Food Additives, Flavourings, Processing Aids and Food Contact Materials (AFC). 122 pp. doi:10.2903/j.efsa.2008.754 

About PD

PD is passionate about applying his background in math, statistics, and economics to apply new and interesting ideas about health, nutrition, and the incentives that drive products and the policies that surround them.