A simple simulation illustrating triplet contacts among different lymphocytes.

simulate_contacts.py

"""Simulates and visualizes triples of DCs, CD4+ and CD8+ T cells."""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.animation as animation
import scipy.spatial as spatial


# Radius around DC within which T cells are considered in contact
contact_radius = 0.1

# (Initial) positions of the cell types
CD8 = np.random.rand(25,2)
CD4 = np.random.rand(10,2)

nDCs = 7
squares = [[x,y]
    for x in range(np.ceil(np.sqrt(nDCs)).astype(int))
    for y in range(np.ceil(np.sqrt(nDCs)).astype(int))
    ]
np.random.shuffle(squares)
DCs = (squares[0:nDCs] + np.random.rand(nDCs,2))/np.sqrt(squares.__len__())

# Set figure w/ axes() for simulation and histograms up
fig = plt.figure(figsize=(8,5))

axes_simulation = plt.axes([0.25/8,0.5/5,4.25/8,4.25/5])
axes_simulation.set_xlabel('Simulation Data')
axes_simulation.axis([0, 1, 0, 1])
axes_simulation.set_xticks([])
axes_simulation.set_yticks([])
CD8_points, = axes_simulation.plot([], [], 'bo')
CD4_points, = axes_simulation.plot([], [], 'ro')

axes_velocity = plt.axes([4.75/8,1/5+2*3.75/3/5,3/8,3.75/3/5])
axes_velocity.set_xlabel('Velocity')
axes_velocity.set_xticks([])
axes_velocity.set_yticks([])
number_measurements = 300
CD8_velocities = []
CD8_velo_hist = [[],[],[]]
CD4_velocities = []
CD4_velo_hist = [[],[],[]]

axes_angle = plt.axes([4.75/8,0.75/5+3.75/3/5,3/8,3.75/3/5])
axes_angle.set_xlabel('Turning Angle')
axes_angle.set_xticks([])
axes_angle.set_yticks([])
CD8_turn_angles = []
CD8_angles_hist = [[],[],[]]
CD4_turn_angles = []
CD4_angles_hist = [[],[],[]]
dx8 = np.random.rand(CD8.shape[0],CD8.shape[1])-1/2
dx4 = np.random.rand(CD4.shape[0],CD4.shape[1])-1/2

axes_contact_times = plt.axes([4.75/8,0.5/5,3/8,3.75/3/5])
axes_contact_times.set_xlabel('Contact Time')
axes_contact_times.set_xticks([])
axes_contact_times.set_yticks([])
existing_contact_times = {}
total_contact_times = []
total_cont_times_hist = [[],[],[]]


def init():
    """Initializes an "empty" frame for FuncAnimation"""
    CD8_points.set_data([], [])
    CD4_points.set_data([], [])
    axes_simulation.plot(DCs[:,0], DCs[:,1], 'yo')
    for DC in DCs:
        axes_simulation.add_patch(plt.Circle(DC, radius=contact_radius,
                                             fc='y', lw=0, alpha=0.33))
    return CD8_points, CD4_points,


def walk(i):
    """Produce changing parts of the plot by random walk"""
    global CD8, CD4
    global CD8_velocities, CD4_velocities, CD8_velo_hist, CD4_velo_hist
    global CD8_turn_angles, CD4_turn_angles
    global CD8_angles_hist, CD4_angles_hist, dx8, dx4
    global existing_contact_times
    global total_contact_times, total_cont_times_hist

    # Silly walk w/ periodic BC according to Nombela-Arrieta et al. 2007
    CD8_previous = np.copy(CD8)
    dx8_previous = np.copy(dx8)
    r = np.random.lognormal(0, 1/2, (CD8.shape[0], 1))/100
    a = (1-np.sign(dx8_previous[:,0])).reshape(-1,1)/2*np.pi + \
        np.sign(dx8_previous[:,0]).reshape(-1,1)*np.arcsin(dx8_previous[:,1]/
        np.sqrt(dx8_previous.dot(dx8_previous.T).diagonal())).reshape(-1,1)
    a += np.random.lognormal(0,1/2, (a.shape[0], 1)) * \
        (2*np.random.randint(0,2,a.shape[0]) - 1).reshape(-1,1)
    CD8 += r*np.concatenate((np.cos(a), np.sin(a)), axis = 1)
    dx8 = CD8 - CD8_previous
    CD8 %= 1

    CD4_previous = np.copy(CD4)
    dx4_previous = np.copy(dx4)
    r = np.random.lognormal(0, 1/2, (CD4.shape[0], 1))/100
    a = (1-np.sign(dx4_previous[:,0])).reshape(-1,1)/2*np.pi + \
        np.sign(dx4_previous[:,0]).reshape(-1,1)*np.arcsin(dx4_previous[:,1]/
        np.sqrt(dx4_previous.dot(dx4_previous.T).diagonal())).reshape(-1,1)
    a += np.random.lognormal(0,1/2, (a.shape[0], 1)) * \
        (2*np.random.randint(0,2,a.shape[0]) - 1).reshape(-1,1)
    CD4 += r*np.concatenate((np.cos(a), np.sin(a)), axis = 1)
    dx4 = CD4 - CD4_previous
    CD4 %= 1

    # Compute velocities & turning angles
    if i < number_measurements:
        CD8_velocities.extend(np.sqrt(dx8.dot(dx8.T).diagonal()),)
        CD4_velocities.extend(np.sqrt(dx4.dot(dx4.T).diagonal()),)

        if i == 1:
            CD8_turn_angles = []
            CD4_turn_angles = []

        dx8_norm = np.sqrt(dx8.dot(dx8.T).diagonal())
        dx8_previous_norm = np.sqrt(
            dx8_previous.dot(dx8_previous.T).diagonal())
        dx4_norm = np.sqrt(dx4.dot(dx4.T).diagonal())
        dx4_previous_norm = np.sqrt(
            dx4_previous.dot(dx4_previous.T).diagonal())

        CD8_turn_angles.extend(
            np.sign(dx8_previous[:,0]*dx8[:,1]-dx8_previous[:,1]*dx8[:,0])*
            np.arccos(dx8.dot(dx8_previous.T).diagonal()/
            dx8_norm/dx8_previous_norm),)
        CD4_turn_angles.extend(
            np.sign(dx4_previous[:,0]*dx4[:,1]-dx4_previous[:,1]*dx4[:,0])*
            np.arccos(dx4.dot(dx4_previous.T).diagonal()/
            dx4_norm/dx4_previous_norm),)

    # Analyze contacts
    CD8_tree = spatial.cKDTree(CD8)
    CD8_contacts = CD8_tree.query_ball_point(DCs, contact_radius)
    CD4_tree = spatial.cKDTree(CD4)
    CD4_contacts = CD4_tree.query_ball_point(DCs, contact_radius)

    current_triples = [
        (iDC, CD8_contact, CD4_contact)
        for iDC, DC in enumerate(DCs)
        for CD8_contact in CD8_contacts[iDC]
        for CD4_contact in CD4_contacts[iDC]
        ]

    for triple in current_triples:
        if triple in existing_contact_times.keys():
            existing_contact_times[triple] += 1
        else:
            existing_contact_times[triple] = 1

    contact_ended = False
    for triple in set(existing_contact_times.keys()) - set(current_triples):
        total_contact_times.append(existing_contact_times[triple])
        contact_ended = True
        del existing_contact_times[triple]

    # Plot cells, contacts and histograms
    triple_triangles = []
    if i > 0:
        for trpl in range(current_triples.__len__()):
            triple_triangles.append(axes_simulation.add_patch(plt.Polygon(
                [DCs[current_triples[trpl][0],:],
                 CD8[current_triples[trpl][1],:],
                 CD4[current_triples[trpl][2],:]],
                fc='g', alpha=0.1)))

    CD8_points.set_data(CD8[:,0], CD8[:,1])
    CD4_points.set_data(CD4[:,0], CD4[:,1])

    if i % 20 == 1 and i < number_measurements:
        CD8_velo_hist = axes_velocity.hist(CD8_velocities, 50, color='blue')
        CD4_velo_hist = axes_velocity.hist(CD4_velocities, 50, color='red')
        CD8_angles_hist = axes_angle.hist(CD8_turn_angles, 50, color='blue')
        CD4_angles_hist = axes_angle.hist(CD4_turn_angles, 50, color='red')

    if contact_ended:
        total_cont_times_hist = axes_contact_times.hist(
            total_contact_times, 50, color='green')

    return tuple(triple_triangles) + (CD8_points,) + (CD4_points,) \
         + tuple(CD8_velo_hist[2]) + tuple(CD4_velo_hist[2]) \
         + tuple(CD8_angles_hist[2]) + tuple(CD4_angles_hist[2]) \
         + tuple(total_cont_times_hist[2])


# To avoid "AttributeError: 'NoneType' object has no attribute 'tk'"
# errors increase interval :>
anim = animation.FuncAnimation(fig, walk, init_func=init, blit=True)

plt.show()

# print(animation.writers.avail)
# mywriter = animation.MencoderFileWriter(fps=10)
# anim.save('simulate_contacts.mp4',writer=mywriter)