Modelling CAR-T cells

Introduction

"How far can a T cell migrate after being activated in a tumor microenvironment?"

Simulation Scenario

CAR-T cells are injected onto the surface of a solid tumor (origin of coordinate system) and are fully activated there.
Through diffusion (modelled as random walk), they leave the tumor microenvironment and enter nearby tissue. T-cells do not attach to any healthy cell, so they keep diffusing.
Once they move away from the tumor microenvironment, they start being deactivated due to lack of inputs that are essential for their survival.

Paremetres

Tumor size (radius): 5 mm
Cells simulated: 4000
Maximum time alive outside microenvironment: 30 hours
Diffusion model: based on speed of T cells, modelled as isotropic, uncorrelated on x, y and z axes (details below)

Fitting Velocity

We model T-cells as moving in "lunges", lasting for 2 minutes each. During each "lunge" they move straight with a speed chosen randomly from a distribution. We fitted a Rayleigh curve to the lateral velocity histogram provided in the paper, as we consider that the cells movements in x, y, and z directions are uncorrelated. Other distributions did not fit the data as good, validating our assumption.

You can find many more details about the modelling in the script iteslf. We believe in Literate programming!

Fitting figure to the histogram of lateral velocity

MatLAB model script


% ETH Zurich
% iGEM 2017 Collaboration
% CAR-T cell model: version 2 [26 Oct 2017]
% Author: Nikolaos Korasidis
%         igem@ethz.ch

%% Cleanup workspace
clear, close all;

%% Model velocity
% Reproduce lateral velocity histogram
% (Updated histogram)
lateral_velocity_counts = [100,75,250,88,432,238,312,238,325,168,212,112,112,88,62,25,25,25,12,12];
interval = 1.5;

% Create synthetic lateral velocity data from the histogram
total_counts = sum(lateral_velocity_counts);
rec_data = zeros(total_counts, 1);
t = 1;
for i = 1:length(lateral_velocity_counts)
counts = round(lateral_velocity_counts(i));
rec_data(t:(t + counts - 1)) = repmat(i * interval - interval/2, counts, 1);
t = t + counts;
end

% Plot histograms and fitted distributions.
figure('NumberTitle', 'off', 'Name', 'Fit Lateral Velocity');
subplot(2, 3, 1);
histfit(rec_data, length(lateral_velocity_counts), 'exponential')
title('Exponential');
subplot(2, 3, 2);
histfit(rec_data, length(lateral_velocity_counts), 'gamma')
title('Gamma');
subplot(2, 3, 3);
histfit(rec_data, length(lateral_velocity_counts), 'lognormal')
title('Lognormal');
subplot(2, 3, 4);
histfit(rec_data, length(lateral_velocity_counts), 'normal')
title('Normal');
subplot(2, 3, 5);
histfit(rec_data, length(lateral_velocity_counts), 'rayleigh')
title('Rayleigh');
subplot(2, 3, 6);
histfit(rec_data, length(lateral_velocity_counts), 'weibull')
title('Weibull');
shg

% Commentary: The best-fitting distributions are Weibull and Rayleigh.  The
% Weibull distribution is actually a generalization of Rayleigh, so we will
% drop it, to prevent overfitting. The Rayleigh distribution arises when we
% take the norm of a vector that has uncorrelated and normally distributed
% gaussian variables as elements. This has a natural interpretation in this
% problem: we are modelling lateral 2D-speed, which arises naturally from
% uncorrelated velocities in the (orthogonal) x and y directions, with the same
% mean (0) and variance.
%
% https://en.wikipedia.org/wiki/Rayleigh_distribution

% Instead of fitting a Rayleigh distribution to the synthetic data, we will use
% the mean velocity value reported in the text, and select a sigma parametre so
% that the resulting Rayleigh distribution has the same mean.
mean_lateral_velocity = 10.7; % um/min
sigmahat = mean_lateral_velocity / sqrt(pi/2);
[mv, v] = raylstat(sigmahat);
fprintf('Mean lateral (x-y) velocity: %.2f um/min and stddev: %.2f um/min\n', mv, v^0.5);

% Given the above, we will assume that, in each lunge, the velocity on the x
% and y axis are independent, and that they are normally distributed random
% variables of mean 0 and standard deviation sigmahat. We will assume that the
% z axis has similar behavior with x and y.

% Thus, to generate a velocity during a lunge, we just sample 3 times from a
% N(0, sigmahat^2) distribution. This automatically gives a velocity vector
% with direction sampled uniformly from the unit ball in R^3 and a
% length(speed) drawn from a hypothetical Rayleigh distribution of modified
% sigmahat, to model spatial, not just lateral speed.
fprintf('Axial velocity along all 3 axes modelled as Gaussian(%.2f, %.2f) um/min\n', 0, sigmahat^2);

% Note: In the paper it is mentioned that the radial velocity (v_z) has a lower
% mean than the lateral velocity (v_x, v_y). We consider this an artifact of
% the experimental procedure, which used very thin tissue layers, thus limiting
% movement in the z axis.

%% Model lunge time.
% In the previous verion, the lunge time was modelled as normally distributed,
% with a mean time 2 min and stddev 0.4 min.
%
% Given that the total distance D travelled per lunge is the product of speed S
% and time T, where S and T are independently drawn, one can show that
%
% Var(ST) = [E(S)]^2 * Var(T) + [E(T)]^2 * Var(S) + Var(S) * Var(T)
%
% Since the variance and the mean of the Rayleigh distribution is much higher
% than that of the normal distribution, we consider the variability in time
% irrelevant and we consider all lunges to last 2mins. This is strongly
% supported by the sharp peak in Figure 3C, suggesting that the modelling as
% normal distribution was too lax to begin with.
%
% https://en.wikipedia.org/wiki/Variance#Sum_of_correlated_variables

% This change also allows us to take time-aligned snapshots of the cells
% modelled and speeds up the simulation considerably.

lunge_duration = 2; % minutes

%% Model time-to-live (TTL) outside of the tumor (OOT)
% In the previous verion, the TTL was modelled as normally distributed, with
% mean 27 hours, and stddev 3 hours. However, this is not worst-case modelling.
% We will shall assume the worst scenario for all cells, that is, a common large
% TTL for all cells.
OOT_TTL = 30 * 60; % minutes

%% Model tumor as spherical
% Smaller tumors make simulation end faster.
tumor_radius = 5 * 1000; % micrometre

%% Simulate
% Setup paremetres.
rng default;
population = 5000;

% For each T-cell, start from a random point at the edge of the tumor.
position = uni3v(population) * tumor_radius;

total_time = zeros(population, 1);
cur_out_time = zeros(population, 1);
net_active_time = zeros(population, 1);
max_distance = zeros(population, 1);

% Go!
tic
parfor i = 1:population
while cur_out_time(i) < OOT_TTL
total_time(i) = total_time(i) + lunge_duration;
% If current position is inside tumor, reset deactivation timer.
if norm(position(i,:)) <= tumor_radius
cur_out_time(i) = 0;
net_active_time(i) = net_active_time(i) + lunge_duration;
else
cur_out_time(i) = cur_out_time(i) + lunge_duration;
end
% Chose a random velocity for the next lunge.
% See modelling section for details.
position(i,:) = position(i,:) + normrnd(0, sigmahat, 1, 3) * lunge_duration;
max_distance(i) = max([norm(position(i,:)) max_distance(i)]);
end
end
toc

%% Plot results
final_distance = sqrt(sum(position.^2, 2));
net_active_fraction = net_active_time./total_time;

total_time_h = total_time / 60;
net_active_time_h = net_active_time / 60;
final_distance_mm = final_distance / 1000;
max_distance_mm = max_distance / 1000;

perc = 0:0.2:1;

figure('NumberTitle', 'off', 'Name', 'Simulation stats');

circle_size = 1.5;

subplot(3, 3, 1);
scatter(total_time_h, final_distance_mm, circle_size);
title('Final distance VS Total Time Alive');
xlim([2e1, 1e4]);
set(gca, 'xscale', 'log'); 
ylabel('[mm]');
ylim([5, 7.5]);
%set(gca, 'yscale', 'log'); 

subplot(3, 3, 2);
scatter(100 * net_active_fraction, final_distance_mm, circle_size);
title('Final distance VS Active Fraction');
ylim([5, 7.5]);
%set(gca, 'yscale', 'log'); 

subplot(3, 3, 3);
histogram(final_distance_mm, 'BinLimits', [5, 7.5], 'BinMethod', 'fd', 'Orientation', 'horizontal');
title('Distance at deactivation');
xlim([0 0.15 * population]);
set(gca, 'XTick', 0.15 * population * perc);
set(gca, 'XTickLabel',sprintf('%1.f\n', 15 * perc));
ylim([5, 7.5]);

subplot(3, 3, 4);
scatter(total_time_h, max_distance_mm, circle_size);
title('Max distance VS Total Time Alive');
xlim([2e1, 1e4]);
set(gca, 'xscale', 'log'); 
ylabel('[mm]');
ylim([5, 7.5]);

subplot(3, 3, 5);
scatter(100 * net_active_fraction, max_distance_mm, circle_size);
title('Max distance VS Active Fraction');
ylim([5, 7.5]);

subplot(3, 3, 6);
histogram(max_distance_mm, 'BinLimits', [5, 7.5], 'BinMethod', 'fd', 'Orientation', 'horizontal');
title('Maximum distance reached');
xlabel('T-cells [%]');
xlim([0 0.15 * population]);
set(gca, 'XTick', 0.15 * population * perc);
set(gca, 'XTickLabel',sprintf('%1.f\n', 15 * perc));
ylim([5, 7.5]);

subplot(3, 3, 7);
histogram(total_time_h, logspace(1, 4, 20), 'Normalization', 'probability');
title('Total Time Alive');
xlabel('[h]');
xlim([2e1, 1e4]);
set(gca, 'xscale', 'log'); 
ylabel('T-cells [%]');
ylim([0 0.25]);
set(gca, 'YTick', perc / 2);
set(gca, 'YTickLabel',sprintf('%1.f\n', 100 * perc / 2));

subplot(3, 3, 8);
histogram(net_active_fraction, 20, 'Normalization', 'probability');
title('Active Fraction');
xlabel('[%]');
set(gca, 'XTick', perc);
set(gca, 'XTickLabel',sprintf('%1.f\n', 100 * perc));
ylim([0 0.25]);
set(gca, 'YTick', perc / 2);
set(gca, 'YTickLabel',sprintf('%1.f\n', 100 * perc / 2));

shg

Results

Here is a summary of the results of the simulation, in grapical format. We can clearly see that the cells reside near the original point of injection for a long time, but they do not move far. Thus, the safety of the CARTEL-TM project is roughly established.

Possible Refinements

To make the model more realistic, one would essentially need to do the following:

Model tumor composition: a tumor should be composed of layers of tissue, with a necrotic zone inside and alive cells on the outer shell.
Model tumor micronvironment: apoptotic tumor cells should release cytokines and chemokines, which diffuse, creating a chemical gradient.
Model atraction and repulsion of CAR-T cells due to chemical gradients
Model attachment of CAR-T cells to tumor tissue

Details concerning T cell migration under the presence of chemical gradients can be found in other resources. [2]

References

Miller, Mark J., et al. "Autonomous T cell trafficking examined in vivo with intravital two-photon microscopy.", Proceedings of the National Academy of Sciences (2003). doi:10.1073/pnas.2628040100
Krummel, Matthew F., Frederic Bartumeus, and Audrey Gérard. "T cell migration, search strategies and mechanisms." Nature Reviews Immunology 16.3 (2016): 193-201.doi:10.1038/nri.2015.16

Team:ETH Zurich/Collaborations/CART-Model