A couple of years ago, Robert had a computer science question. One at the data analysis and visualization stage of the relationship between computer science and cancer. Given that I haven’t posted code on TheEGG in a long time, I thought I’d share some visualizations I wrote to address Robert’s question.

There are many ways to measure the size of populations in biology. Given that we use it in our game assay, I’ve written a lot about using time-lapse microscopy of evolving populations. But this isn’t the only — or most popular — approach. It is much more common to dillute populations heavily and then count colony forming units (CFUs). I’ve discussed this briefly in the context of measuring stag-hunting bacteria.

But you can also combine both approaches. And do time-lapse microscopy of the colonies as they form.

A couple of years ago, Robert Vander Velde Andriy Marusyk were working on experiments that use colony forming units (CFUs) as a measure of populations. However, they wanted to dig deeper into the heterogeneous dynamics of CFUs by tracking the formation process through time-lapsed microscopy. Robert asked me if I could help out with a bit of the computer vision, so I wrote a Python script for them to identify and track individual colonies through time. I thought that the code might be useful to others — or me in the future — so I wanted to write a quick post explaining my approach.

This post ended up trapped in the drafts box of TheEGG for a while, but I thought now is as good a time as any to share it. I don’t know where Robert’s work on this has gone since, or if the space-time visualizations I developed were of any use. Maybe he can fill us in in the comments or with a new guest post.

So let’s just get started with the code.

Of course, we first need to import the main packages: numpy, pyplot, and cv2.

import numpy as np
import cv2
from matplotlib import pyplot as plt

The first two are standard packages, the last one — OpenCV — takes a little bit more work to install.

Now we do two main tasks at once, we load all the images and create something I want to call a ‘space-time map’. A space-time map is an image that uses the colour map of a pixel to represent the number of time points that it appears in. This is the first name that occurred to me, if you’ve seen this visualisation used before, dear reader, and know its name then please let me know.

threshImgs_all = []
num_imgs = 24

#load the images and create space-time image (total_img)
img_all = []
total_img = np.zeros_like(cv2.imread(str(0).zfill(4) + '.tif',cv2.IMREAD_GRAYSCALE))

for img_num in range(0,24):
    f_name = str(img_num).zfill(4) + '.tif'
    img = cv2.bitwise_not(cv2.imread(f_name,cv2.IMREAD_GRAYSCALE))
    img_all.append(img)
    total_img = cv2.scaleAdd(img,1/num_imgs,total_img)

plt.imshow(total_img, cmap="magma_r")
plt.show()

This results in an image like:

From this image, we can get persistent numbers for all the colonies that existed:

#get the colonies
_, total_thresh = cv2.threshold(img,127,255,cv2.THRESH_BINARY)
_, total_colonies = cv2.connectedComponents(total_thresh)

num_colonies = np.amax(total_colonies)
print("There are " + str(num_colonies) + " colonies")

More importantly, the image total_colonies now has each non-background pixel labeled by its colony number, so counting the number of pixels in each colony at each time point becomes as straightforward as applying a mask:

#use the total image (properly thresholded) as the permanent numbers for the colonies; 
#get future colonies numbers from them)

colony_sizes = np.zeros((num_colonies + 1,num_imgs), dtype=np.int) 
#Note that colony_size[0,:] will contain the amount of empty space

img_num = 0
for img in img_all:
    #label colonies by their numbers (for upto 255 colonies):
    labeled_img = np.minimum(total_colonies,img)
    
    #get the colonies that appear and their sizes
    colonies, sizes = np.unique(labeled_img, return_counts = True)
    colony_sizes[colonies,img_num] = sizes
    img_num += 1

#plt.imshow(total_colonies, cmap='magma_r')
for colony in range(1,num_colonies + 1):
    plt.plot(colony_sizes[colony,:])

plt.yscale('log')
plt.show()

Unfortunately, there is a number of colonies that ‘blink’ in and out of existance. This is not a manifestation of reality, but probably an artefact of the image processing software used to produce the initial threshold images and the sensitivity of the microscope. As such, it can be helpful to clean up the time series and focus only on the colonies that didn’t go to extinct during the experiment and look at their population dynamics.

#let's clean up by eliminating colonies that go extinct at some point.

colony_lifetimes = np.sum(colony_sizes > 0, axis = 1)
surviving_colonies = np.where(colony_lifetimes == num_imgs)[0][1:]

print(surviving_colonies)

for colony in surviving_colonies:
    plt.plot(colony_sizes[colony,:])

plt.yscale('log')
plt.show()

But the figure that this produces still a difficult figure to make sense. I don’t even bother to produce it given how many different lines there are going in different directions for growth rates.

What we really care about is higher level properties of these colonies like their growth rate, so let’s infer those with the help of scipy:

#among those that don't go extinct, let's calculate the growth rates

from scipy.stats import mstats

growth_rates = np.zeros(num_colonies + 1)
growth_rates_low = np.zeros(num_colonies + 1)
growth_rates_high = np.zeros(num_colonies + 1)

for colony in surviving_colonies:
    growth_rates[colony], _, growth_rates_low[colony], growth_rates_high[colony] = \
        mstats.theilslopes(np.log(colony_sizes[colony,:]))

plt.errorbar(np.arange(num_colonies + 1),growth_rates,
             yerr = [growth_rates - growth_rates_low, growth_rates_high - growth_rates],fmt='ko')
plt.xlabel('Colony number')
plt.ylabel('Growth rate')

plt.show()

This yields an easier to look at colony growth rate plot, with 95% confidence intervals.

Above, we have a fitness measure for each colony, so we can look not only at the number of colony forming units but also at differences in how the colonies formed. I still find it hard to make sense of this particular plot, but looking explicitly at the inter-colony heterogeneity does seem like a good exercise. Definitely better than just summarising it as a single variance. Especially since I know from experience that sometimes a variance can hide an interesting discovery.

How would you, dear reader, extend these visualizations? Or is there a good use that you can think of putting them to? After all, visualizations are one of the most important parts of science. I hope this code helps a little. At least as inspiration, or an example of how easy it is to get things done with Python.

The motivation for this post came about a month ago when Nathan Farrokhian was asking for some advice on how to repeat our game assay with a new experimental system. He has since done so (I think) by measuring the game between Gefitinib-sensitive and Gefitinib-resistant cell types. And I thought it would make a nice post in the quick introductions series.

Of course, the details of the system don’t matter. As long as you have an array of growth rates (call them yR and yG with corresponding errors yR_e and yG_e) and initial proportions of cell types (call them xR and xG) then you could repeat the assay. To see how to get to this array from more primitive measurements, see my old post on population dynamics from time-lapse microscopy. It also has Python code for your enjoyment.

In this post, I’ll go through the two final steps of the game assay. First, I’ll show how to fit and visualize fitness functions (Figure 3 in Kaznatcheev et al., 2019). Second, I’ll transform those fitness functions into game points and plot (Figure 4b in Kaznatcheev et al., 2019). I’ll save discussions of the non-linear game assay (see Appendix F in Kaznatcheev et al., 2019) for a future post.

As always, it is good to start with some imports and visualize our data:

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import mstats

r_GC = 'c' #colour for R-type cells
g_GC = 'm' #colour for G-type cells

plt.errorbar(xR,yR,yerr=yR_e,fmt=r_GC+'s',markersize=ms,fillstyle='none')
plt.errorbar(xG,yG,yerr=yG_e,fmt=g_GC+'s',markersize=ms,fillstyle='none')

At this point, you’ll have a plot of growth-rates versus initial proportion. This is what you need to get fitness functions. And at this point your first choice enters into the analysis. Based on your knowledge of the system and the goal of your measurement, you should decide on what sort of class of functions should be considered as candidates for the fitness functions. I’ll go into details of this in the future post, but if our goal is to measure matrix (i.e. linear) games then we should choose linear functions.

So, let’s fit some linear functions and visualize.

[ra,rb], [[rav,rabv],[_,rbv]] = np.polyfit(xR,yR,1,w=1/yR_e**2,full=False,cov=True)
[ga,gb], [[gav,gabv],[_,gbv]] = np.polyfit(xG,yG,1,w=1/yG_e**2,full=False,cov=True)

plt.plot([0,1],[rb,rb + ra],r_GC+'--')
plt.plot([0,1],[gb,gb + ga],g_GC+'--')

If you were to do this with our data for Alectinib-sensitive versus Alectinib-resistant non-small-cell lung cancer. And repeated it four times for four conditions (DMSO, DMSO + Fibroblasts, Alectinib, Alectinib + Fibroblasts) — plotting each condition in a new subfigure then (with a few extra axis labels, etc) you’d arrive at a graphic like the following:

Now, the last step is to visualize these fitness functions as points in game-space. For this, we want to propagate our errors, so it’s helpful to define the following auxiliary functions:

def FitRightPoint(fit):
out = fit[0] + fit[1]
out_v = fit[2] + fit[3] + 2*fit[4]

return [out, out_v]

def IndepDiff(point1,point2):
out = point1[0] - point2[0]
out_v = point1[1] + point2[1]

return [out, out_v]

def GameCoord(fit_L,fit_R):
point_y = IndepDiff(fit_R[[1,3]],fit_L[[1,3]])
point_x = IndepDiff(FitRightPoint(fit_L),FitRightPoint(fit_R))

return point_x, point_y

Now, you can just calculate the game-point and plot:

[game_x,gxv],[game_y,gyv] = ...
GameCoord(np.array([ra, rb, rav, rbv, rabv]), np.array([ga, gb, gav, gbv,gabv]))
game_xe = np.sqrt(gxv) #convert variance to error
game_ye = np.sqrt(gyv) #convert variance to error

plt.figure() #make a new figure, to avoid overwriting our fitness function plot

plt.errorbar(game_x,game_y,xerr = game_xe, yerr = game_ye, fmt='k.')

Again, if you repeat this four times and add some cosmetic colours for the points and labels then you’ll end up with a nice game space plot. At right, you can see what that looks like for our non-small-cell lung cancer system.

I think that any measurement of a game is interesting — since we have so few in cancer. However, the measurements are even more exciting if game points for different conditions end up in different quadrants on the game space. Each quadrant corresponds to a qualitatively different kind of dynamic, so seeing games in different quadrants means that we’re seeing a switch in the type of interaction. This can be especially important if we’re aiming to treat the game instead of the player.

The most interesting quadrants are probably the top right and bottom left, where internal fixed points exist. The fixed point is fixed in the top right and unstable in the bottom left. I’ve yet to measure a system with a game in the bottom-left corner. I’ll be sure to report back when we find out.

If you measure some evolutionary games with the assay, dear reader, then please let me know!

And I hope this code sketch is useful. I’ll soon post more detailed game assay code on GitHub.

References

Kaznatcheev, A., Peacock, J., Basanta, D., Marusyk, A., & Scott, J. G. (2019). Fibroblasts and Alectinib switch the evolutionary games played by non-small cell lung cancer. Nature Ecology & Evolution, bioRxiv: 179259.

In this case, I opted to remove all mentions of the bone and cancer. Instead, I asked a simpler but more abstract question: what qualitative features must a minimal model of the evolution of resistance have in order for drug holidays to be superior to a single treatment block? In this post, I want to set up this question precisely, show why drug holidays are difficult in evolutionary models, and propose a feature that makes drug holidays viable. If you find this topic exciting then you should consider registering for the 6th annual Integrated Mathematical Oncology workshop at the Moffitt Cancer Center.^[1] This year’s theme is drug resistance.

Discontinuous treatment powered by holidays

Most discussions of treatment holidays are focused on managing toxicity. Especially in the case of chemotherapy, the treatment is often devastating to the patient. For many therapies, the limiting factor is for how long (or how strong) a patient can be treated without being killed by the therapy. Thus, these holidays aren’t the active ingredient of the therapy, they are not what lets the therapy supress cancer cells. These holidays are a way to manage side-effects. And although I think it can be very important to focus on the toxicity of treatment, in this case I want to think about the case where holidays are a central driver behind the effectiveness of therapy. Is there a class of therapies and evolutionary dynamics such that the drug fails to treat the tumour unless its application is interspersed with holidays of no treatment?

In particular, I want to consider two therapies of the same strength (per unit time) that have the same start and end time. The continuous case would run continuous treatment from the start time to the end time. The discontinuous therapy punctuates blocks of treatments with holidays from the start time until the end time. Thus, I want to consider cases where the discontinuous therapy has strictly less therapy. And yet, I want to find settings under which the tumour burden at the end time is lower for the discontinuous therapy.

Minimal model and the difficulty of discontinuous treatment

Of course, in a real clinical setting, we might not care if a discontinuous treatment works by allowing a different dosage, recovery of sensitivity, or booting up the patient’s immune system. We only really care that it works. But if we consider all three (or more) effects at once in a heuristic model — or in a big simulation — then we won’t really understand why our strategy succeeds. Thus, in this post, I want to isolate just the recovery-of-sensitivity-during-treatment-holidays aspect of the strategy.

I will start with one of the simplest possible model of resistance. A tumour (or pathogen) growing towards carrying capacity. There is a single therapy and two tumour subtypes: a sensitive type of density that has one fitness in the absence of treatment and a second fitness in the presence of treatment; and a resistant type of density that a fitness regardless of if therapy is present of absent. Giving us the dynamics in the presence of therapy:

where is the total tumour burden, and where we remove the hat over for the times that treatment is off. On the one hand, it is easy to see that we have no hope at all if or if . On the other hand, if then the tumour is cured with continuous therapy. Thus, we will focus on .

Unfortunately, under these conditions, it is clear that for all time whenever and . Thus, the density of resistant cells will always be increasing and will only slow down as total tumour burden approaches 1.

Resensitization of the resistant population

To overcome this, we need a way for . This can be done by allowing for a flow from resistant to sensitive cells in the absence of therapy. In other words, when therapy is absence, our equations will be:

There are several ways we could interpret these equations, but my favourite option is to think of them as individual cells losing resistance. When under stress, each resistant cell re-evaluates its commitment to investing energy in resistance and gives up that strategy at a rate of . It is important to note, that unlike mutation, this is not linked to reproduction. In particular, does not necessarily need to depend on the difference , and — more importantly — it is independent of the remaining free space .^[2] This means that we can have the suppression of resistance when therapy is off and . In other words, there is hope of controlling the resistant population in large tumours, even if we cannot eliminate the whole population.

When treatment is on, the resistant proportion will grow back but it is possible to still have the total tumour burden shrink by killing sensitive cells quickly enough. In order for total tumour volume to shrink, we need which is achieved only when with where is the fraction of the tumour that is resistant. Thus, the tumour will shrink whenever the fraction of sensitive cells is high enough; i.e. when . Of course, this means that the fraction of sensitive cells will also decrease and once it falls below the threshold, we need to repeat the cycle of turning therapy on and off.

Unfortunately, just those conditions are not enough to guarantee that discontinuous therapies can control the tumour burden and outperform a single continuous treatment. The rates of loss of sensitivity under treatment and growth of tumour burden without treatment have to be low enough to be surpassed by the rate of loss of resistant density without treatment and loss of tumour density with treatment. It is probably possible to calculate these rates analytically and also solve for the average burden of the controlled tumour. However, this was meant to be a quick post so I only present a proof of concept by simulation.³

Proof of concept and adaptive therapy

Below is a figure showing an example of when discontinuous therapy outperforms a continuous therapy of equal strength. All parameters — except the times that treatment is turned on or off — are the same across the two simulations of the ODEs. These parameters are: and the initial conditions . This gives a critical proportion of sensitive cells , above this tumours shrink under treatment and below it they grow. This critical proportion is shown in solid green; in dotted green is the actual proportion of sensitive cells . In solid black is the total tumour burden and in red is the resistant tumour burden . For the times that therapy is off the background is white, and when therapy is of the background is gray.

In the top panel with a continuous therapy, as expected, the resistant population quickly takes over the whole tumour, grows the tumour to carrying capacity, and leaves the patient with a tumour burden of 1. With discontinuous treatment, it is possible to time treatment holidays such that the resistant tumour burden is kept oscillating around ~0.3 and the total tumour burden is kept oscillating around ~0.85 without ever surpassing ~0.93. This is not a huge reduction in tumour burden but it does show how holidays can help us control a tumour that is unresponsive to continuous therapy.

If you are curious, the discontinuous treatments in the second panel are from times 1 to 2, 3.5 to 4.5, 6 to 9, 10 ti 13, 14 to 17.5, 18.5 to 22.5, 23.5 to 27, 28 to 32, and 33 to 35. While the continuous treatment is just a single block from 1 to 35. In other words, the discontinuous treatment includes a total of 9 time-steps of holidays distributed over eight holidays (the first two are 1.5 time-steps long while the rest are 1). Thus, the discontinuous treatment produces a better outcome even though it ends up applying less than 74% of the drug-hours of the continuous treatment.

In this example, I tuned the treatment timings by eye but it is not hard to adapt my strategy into an adaptive therapy. A theoretical route might be to pick up upper sensitivity threshold . Whenever reaches (or exceeds) turn on the treatment and when reaches (or falls below) turn off the treatment. In practice, though, distinguishing sensitive from resistant cells to measure the proportion might be difficult. Instead, we could just track the tumour burden. We only care about crossing because it changes the sign of under treatment. For a more practical route, we might turn on treatment whenever the tumour burden passes a threshold , and if treatment is on we will turn off the treatment if the tumour burden stopped decreasing.

The practical route is a pretty obvious strategy for adaptive therapy. And I wouldn’t be surprised if doctors do something like this already. Although usually, they switch to a new drug when a tumour stops responding to an existing one. Thus, the difficulty becomes characterising the sort of tumours (or evolutionary dynamics) for which this discontinuous strategy will outperform a continuous counterpart. Then we can know when to return to and old drug to achieve tumour control instead of switching to a new drug in the hopes of a cure.

Notes

1. Whether you work in mathematical oncology or in more general modeling of biological systems then I recommend this event. I have attended for the past three years and plan to attend this year. I greatly enjoyed each prior year: 2013, 2014, 2015. Here is a video of the 3rd workshop to whet your appetite:
2. Note that the structural feature of fitness-based growth rate depending on but the phenotype-switching being independent of it is essential for the effect observed in this model. That means that logistic growth is not an accidental but an essential part of this model. If we want the same sort of results with exponentially (instead of logistically) growing populations then it will have to depend on a different sort of mechanism.
3. I would prefer to not succumb to the curse of computing and avoid relying on ad hoc examples of parameters from simulations. From playing with simulations, it looks like the result is robust. However, for the parameters I tried, the average burden of the controlled tumour is still uncomfortably high. It seems that each parameter setting has an associated region of control, but I cannot figure out how tweaks to the parameters exactly affect the average controlled burden. An analytic treatment would help me resolve this. Unfortunately, I can only come up with rough analytic approximations by discretizing the cycles into two steps. This is not unreasonable — especially since under log or logit transform the ‘waves’ look close to piecewise linear — but not perfect. If you have ideas on a more grounded analytic treatment, dear reader, then let me know. I will also play around on my own when I have time and report back with any successes.

Above, you can see two types of non-small cell lung cancer cells growing in the presence of 512 nmol of Alectinib. The cells fluorescing green are parental cells that are susceptible to the drug, and the ones in red have an evolved resistance. In the 3 days of the video, you can see the cells growing and expanding. It is the size of these populations that we want to quantify.

In this post, I will remedy last week’s omission and share some empirical population dynamics. As before, I will include some of the Python code I built for these purposes. This time the code is specific to how our microscope exports its data, and so probably not as generalizable as one might want. But hopefully it will still give you some ideas on how to code analysis for your own experiments, dear reader. As always, the code is on my github.

Although the opening video considers two types of cancer cells competing, for the rest of the post I will consider last week’s system: coculturing Alectinib-sensitive (parental) non-small cell lung cancer and fibroblasts in varying concentrations of Alectinib. Finally, this will be another tools post so the only conclusions are of interest as sanity checks. Next week I will move on to more interesting observations using this sort of pipeline.

Batch processing the images

Like last time, we are starting with raw 3-channel image data that we need to digitize and convert into the size of populations via our AreaCount method. Our microscope exports these images as an archive directory of files that we need to navigate to process them one at a time. This is where Python shines as a scripting language with easy file operations and directory navigation via the os library.

In particular, we will follow a depth-first traversal of our ScanData directory, looking for the subdirectories that we know to house the images from each time step. After processing the images, we will save a processed copy of the image for later inspection and use, and save the populations to an array for plotting and other further analysis:

import FluorescentArea as fa
import os
import time
import numpy as np

def ConvertArchive(root_dir, subdir_match, col_let_array, row_num_array):
    red_list_all = []
    green_list_all = []
    
    for dirName, subdirList, fileList in os.walk(rootDir):
        if len(subdirList) > 0:
            if subdirList[0] == subDirMatch:
                #track time
                start_time = time.clock()

                #create slices for this timestep
                red_slice = np.zeros((len(col_let_array),len(row_num_array)))
                green_slice = np.zeros((len(col_let_array),len(row_num_array)))

                for col_let_num in range(len(col_let_array)):
                    col_let = col_let_array[col_let_num]
                    for row_num_num in range(len(row_num_array)):
                        row_num = row_num_array[row_num_num]
						
						#analyze each of the three images
                        nums_1, im = fa.AreaCount(col_let,row_num,1,
							dirName+'/'+subDirMatch+'/')
                        cv2.imwrite(dirName + '/'+subDirMatch+ \
							'/yol_'+col_let+str(row_num)+'-1.jpg',im)

                        nums_2, im = fa.AreaCount(col_let,row_num,2,
							dirName+'/'+subDirMatch+'/')
                        cv2.imwrite(dirName + '/'+ subDirMatch+ \
							'/yol_'+col_let+str(row_num)+'-2.jpg',im)

                        nums_3, im = fa.AreaCount(col_let,row_num,3,
							dirName+'/'+subDirMatch+'/')
                        cv2.imwrite(dirName + '/'+subDirMatch+ \
							'/yol_'+col_let+str(row_num)+'-3.jpg',im)

                        red_slice[col_let_num,row_num_num] \
							= nums_1[0] + nums_2[0] + nums_3[0]
                        green_slice[col_let_num,row_num_num] \
							= nums_1[1] + nums_2[1] + nums_3[1]

                red_list_all.append(red_slice)
                green_list_all.append(green_slice)

                end_time = time.clock()

                print(dirName + ' took ' \
					+ str(end_time - start_time) + ' seconds')
	
	return red_list_all, green_list_all

This can then easily be run on any of our 3-channel, 3-positions per well microscopy archives. For example, for Jeff’s first 96-well plate, I would run:

red_all, green_all = ConvertArchive('./ScanData', '325', 
    ['A','B','C','D','E','F'], range(13))

#write the output to file for later use
np.save('./ScanData/red_area',red_all)
np.save('./ScanData/green_area',green_all)

Plotting and fitting growth rates

With all the data recorded, we can start to plot and estimate growth curves. The whole experiment is kept in exponential phase, avoiding confluence, so we will stick to semilog plots. In these, exponential growth becomes straight lines, and thus very easy to measure. For a change of pace, we will actually consider fitting a piecewise linear model with two pieces to semilog data (although the code can fit single lines with the flag flip pwl = False).

import matplotlib
import matplotlib.pyplot as plt
from scipy import optimize
	
def AlectinibPlot(r_1,g_1, y_len, step_size = 1, pwl = True ):
    matplotlib.rcParams['figure.figsize'] = (20, 10)
	
	if pwl: #if we are fitting piecewise linear then define the model
		def piecewise_linear(x, x0, y0, k1, k2):
			return np.piecewise(x, [x < x0], 
				[lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])
		red_rates = []
		green_rates = []
	else: #if fitting linear rates, then initialize output array
		red_rates = np.zeros((y_len,2))
		green_rates = np.zeros((y_len,2))
	
	#initialize our arrays of plot axes
    f, plot_arr = plt.subplots(2,y_len,sharex='all', sharey='all')
    plt.yscale('log')
	
    for y in range(y_len):
        for offset in range(2):
			
			#create all_t to collect data from all 3 replicates
            x_all_t = []
            r_all_t = []
            g_all_t = []
			
            for i in range(3):
                x_all_t = np.append(x_all_t,
					range(len(r_1[::step_size,offset*3 + i,y])))
                r_all_t = np.append(r_all_t,r_1[::step_size,offset*3 + i,y])
                g_all_t = np.append(g_all_t,g_1[::step_size,offset*3 + i,y])
				
				#plot the red and green dynamics
                plot_arr[offset,y].plot(range(len(r_1[::step_size,offset*3 + i,y])),
					r_1[::step_size,offset*3 + i,y],'r-',
                    range(len(g_1[::step_size,offset*3 + i,y])),
					g_1[::step_size,offset*3  + i,y],'g-')

            if pwl: #fit piecewise linear of two lines
                r_p , r_e = optimize.curve_fit(piecewise_linear, x_all_t, 
					np.log(r_all_t))
                g_p , g_e = optimize.curve_fit(piecewise_linear, x_all_t, 
					np.log(g_all_t))
				
				#plot the lines of best fit
                xd = np.linspace(0, len(r_1[::step_size,0,y]), num=100)
                plot_arr[offset,y].plot(xd,np.exp(piecewise_linear(xd,*r_p)),'k-', 
					xd,np.exp(piecewise_linear(xd,*g_p)),'k-')
            else: #calculate the single line of best fit:
                [ra,rb] = np.polyfit(x_all_t,np.log(r_all_t),1)
                [ga,gb] = np.polyfit(x_all_t,np.log(g_all_t),1)
				
				#record the growth rates
				red_rates[y,offset] = ra
				green_rates[y,offset] = rb
				
				#figure out the line starting and ending points
                r_start = np.exp(rb)
                r_stop = np.exp(ra*len(r_1[::step_size,0,y]) + rb)

                g_start = np.exp(gb)
                g_stop = np.exp(ga*len(g_1[::step_size,0,y]) + gb)
                
				#plot the lines from their start/end points
                plot_arr[offset,y].plot([0,len(r_1[:,0,y])],[r_start,r_stop],'k-', 
					[0,len(g_1[:,0,y])],[g_start,g_stop],'k-')
	
	return red_rates, green_rates

So we can run it on our data from before — AlectinibPlot(red_all,green_all,12) — to produce the following graph (with some extra cosmetics added in paint):

Dynamics of fluorescent area of alectinib-sensitive non-small cell lung cancer cells (parental; red) and fibroblasts (green) in different concentrations of Alectinib. Each panel is plotted on a semilog axis, with the x-axis as time running over approximately 6 days, and contains data from 3 separate plates — one for each line. The line for each plate is made from processing 3 images at different positions in the plate. In black is the best fit possible with a function of two piecewise lines. Alectinib concentration varies by panel from left to right: 0, 0.5 nanomols, 1 nM, and then doubling in each successive panel until 512 nM. The top and bottom panels are identical, except in their initial proportion of fibroblasts, with the top having a significantly higher initial proportion.

Alternatively, we can run the above with linear fits by setting the flag pwl = False, and plot the resultant rates versus concentration of alectinib (see below-right). From these two figures we can conclude some obvious things as part of a sanity-check of my analysis:

Growth rate of parental cells (red) and fibroblasts (green) vs. concentration of alectinib. Crosses for high initial number of fibroblasts (top panels of previous figure) and dots to low numbers (bottom panels).

• For most of the conditions (38/48 lines of best fit), a constant exponential growth rate is a reasonable assumption. This suggests that we can use the growth rate parameter as a measure of fitness to feed into later models.
• The parental cells are indeed sensitive to Alectinib, seeing a drastically lower — although still positive — growth rate at higher concentration of the drug.
• The final size of the fibroblast population is not strongly affected by the drug (completely unaffected at concentrations lower than 64 nM). When it is affected, it seems that this is more likely mediated by numbers of fibroblasts and parental cells. This suggests that the fibroblasts need some threshold of signal from parental cells per fibroblast to trigger growth.
• The growth rate of parental cells is not significantly affected by low/high initial proportion of fibroblasts. This suggests that the public good that we expect fibroblasts to be producing — hepatocyte growth factor — is non-linear; maybe a threshold good.

From other evidence, however, we know that the growth rates of the parental cells is much higher in the presence of fibroblasts. Without them the growth of the alectinib-sensitive cells would have been negative, instead of just slightly positive, in the higher drug concentrations. However, this cannot be read from the data in this post. The goal this week is just to spell out some more tools, so that next week I can share some of the accidental science that I presented yesterday at the 2016 Moffitt Scientific Symposium.

To move further into evolutionary game theory, we will need to analyze a different set of experiments. The ones discussed here are great for getting a grip on the biology, and calibrating the analysis pipeline, but there are simply too many games going on. Each of the 12 concentrations of Alectinib are effectively a different game, giving us only 2 to 4 data-points per game. Next week, we will limit ourselves to just two concentrations of Alectinib (0 vs. 512 nM), and expand to all three cell types: alectinib-sensitive, resistant, and fibroblasts.

Together with David Basanta and Jacob Scott, the five of us are looking at the in vitro dynamics of resistance to Alectinib in non-small cell lung cancer. Alectinib is a new ALK-inhibitor developed by the Chugai Pharmaceutical Co. that was approved for clinical use in Japan in 2014, and in the USA at the end of 2015. Currently, it is intended for tough lung cancer cases that have failed to respond to crizotinib. Although we are primarily interested in how alectinib resistance develops and unfolds, we realize the importance of the tumour’s microenvironment, so one of our first goals — and the focus here — is to see how the Alectinib sensitive cancer cells interact with healthy fibroblasts. Since I’ve been wanting to learn basic computer vision skills and refresh my long lapsed Python knowledge, I decided to hack together some cell counting algorithms to analyze our microscopy data.^[2]

In this post, I want to discuss some of our preliminary work although due to length constraints there won’t be any results of interest to clinical oncologist in this entry. Instead, I will introduce automated microscopy to computer science readers, so that they know another domain where their programming skills can come in useful; and discuss some basic computer vision so that non-computational biologists know how (some of) their cell counters (might) work on the inside. Thus, the post will be methods heavy and part tutorial, part background, with a tiny sprinkle of experimental images.^[3] I am also eager for some feedback and tips from readers that are more familiar than I am with these methods. So, dear reader, leave your insights in the comments.

Automated microscopy

A few years ago, I wrote about droids taking academic jobs from a typical theorist’s focus on automating literature review and writing. However, the real automation breakthroughs are on the experimental side. There used to be a time when people would have to look through microscopes by eye, count cells manually, and sketch what they see. I think I even did that at some point in school. Now, Jeff can prepare a 96-well plate and put it in a machine that will take precise photos of the same exact positions in each well, every 2 hours.^[4] Six days later, I can ssh into the same machine and get about 51 thousand images — 40 GB of TIFs — to analyze.

However, even under hundredfold magnification, cells are hard to see and tell apart. So a number of physical and biological tricks are used before the images are even recorded. The first of these is phase-contrast. In my grade-school experience with microscopes, I had a light source under the sample, it would shine through it, with some of the light being absorbed (or otherwise obstructed) by the cell, and thus creating darker patches for me to identify as life. This is called bright-field microscopy, and only takes advantage of the amplitude of the light — that which makes things brighter or darker to our eyes and cameras. But light carries a second channel: phase. Phase is shifted by passing through the different media of the cell membrane and other components, with the amount of shift depending on thickness and refractive index of the object. This information is then captured by the digital camera, combined with the standard bright field information, and converted into an optical thickness for each pixel. The result is a phase-contrast image that is much clearer than standard bright-field and that allows us to see details that would normally only be visible through staining of dead cells.

But even if you can see the cells, it is still often impossible to tell them apart because different cell types can be morphologically identical. To overcome this, a biological engineering trick can be combined with fluorescence imaging. We can take our separate cell lines — in our case: parental cancer cells, resistant cancer cells, and fibroblasts — and add an extra piece of DNA to code for a fluorescing (but otherwise non-functional) protein. Then, during image capture, we shine a light of the excitatory wavelength that is absorbed by the protein and then re-emitted at a longer wavelength which is recorded (after the source wavelength is physically filtered out) by the digital camera. Thus, instead of looking at obstructions (and phase-shift) of an external light source, we end up using the protein itself as the light source and recording where that light is brightest. By giving different cell lines different proteins that are excited by and emit different wavelengths, we can then tell those cells apart by the wavelength of their fluorescence.

The practical result is that for each snapshot of part of a well,^[5] we end up with three photos: one phase-contrast, and one color channel for each of the two fluorescent proteins. The programming problem becomes transforming these 3 photos into a measurement of size for the populations of cancer cells and fibroblasts.

An example of the 3 channels that are recorded from microscopy. From left to right: phase contrast, GFP, mCherry. All data is intensity, but I used the green and red colour space for the 2nd and 3rd image to make them easier to distinguish at a glance.

In our experimental system, we created two variants of each of our cancer and fibroblast cell lines by inserting the genes coding for GFP (absorption 395 nm, emission 509 nm) into one and mCherry (587/610 nm) into the other. These genes are inserted throughout the genome, and so we expect them to be expressed and thus to fluoresce whenever proteins are being transcribed. Further, our fluoroproteins are not confined to any specific part of the cell and distribute themselves haphazardly throughout the cell. This allows the fluorescent area to serve as a good unit of size for our populations.

Conveniently — and not at all by coincidence — fluorescent area is much easier to measure than the number of cells. Mostly, this is due to the tendency of the cancer cells to clump, making segmentation difficult, and the lanky and faint fibroblasts making it hard to find their boundaries.^[6] Of course, this means that as we eventually transform this into replicator dynamics, we will have be mindful of our units being fluorescent area and not cell count, but that choice of units makes no difference to the mathematics.

Basic image processing

“If you want to make an apple pie from scratch,” said Carl Sagan in Cosmos (1980), “you must first invent the universe.” This is equally true for programming. If I want to build my analysis code ‘from scratch’ then I can go arbitrarily far back to progressively more basic building blocks. At some point, if I don’t want to ‘invent the universe’, I have to pick a set of packages to start from. I decided to start from the Python equivalent of store-bought crust and canned apple pie filling: the OpenCV package (and numpy, of course).^[7]

import numpy as np
import cv2

OpenCV makes file operations with images extremely simple:

def AreaCount(col_let, row_num, fot_num, dirName = '', ret_img = True):
    #load the images
    head = dirName + col_let + str(row_num) + '-' + str(fot_num)
    
    rdImage = cv2.imread(head + '-C1.tif',cv2.IMREAD_UNCHANGED)
    gnImage = cv2.imread(head + '-C2.tif',cv2.IMREAD_UNCHANGED)
    if ret_img:
        grImage = cv2.imread(head + '-P.tif',cv2.IMREAD_UNCHANGED)

In this case, the outputs of the cv2.imread functions are simply two-dimensional arrays of intensities. Although the microscope itself creates different kinds of TIFs for the phase-contrast and fluorescent photos, so the output arrays require a conversion into a common 8-bit format:^[8]

    #switch to 8bit with proper normalizing
    if np.amax(rdImage) > 255:
        rdImage8 = cv2.convertScaleAbs(rdImage, alpha = (255.0/np.amax(rdImage)))
    else:
        rdImage8 = cv2.convertScaleAbs(rdImage)
        
    if np.amax(gnImage) > 255:    
        gnImage8 = cv2.convertScaleAbs(gnImage, alpha = (255.0/np.amax(gnImage)))
    else:
        gnImage8 = cv2.convertScaleAbs(gnImage)

To go beyond file operations, thought, I need a brief detour into the real apples of OpenCV.

Straight out of the box, OpenCV gives me tools like histogram equalization for enhancing image contrast. Histogram equalization for contrast works on the premise that what makes an image high contrast is a maximization of information conveyed by the intensity of each pixel. In other words, we want to stretch out the intensity histogram into the highest-entropy version: a uniform distribution. This is a relatively straightforward process, just take the cumulative distribution for the pixels, and then create a function such that a pixel of value with maps to a new value , that way in the new image . Not a hard function to implement for yourself, but with OpenCV it is just one line: img_hicon = cv2.equalizeHist(img).

An example of histogram equalization that makes obvious the colour gradient present in the mCherry channel. On the right of each image is the histogram for their intensity distribution, with the raw histogram in red and the cumulative distribution as a black line. For the equalized image, I am using a different color map with blue corresponding to low intensity and red to high intensity. In the bottom histogram, I am also including bigger bins to make it clearer how the distribution is moved towards a uniform one.

What you might notice right away is the colour gradient radiating outwards. This is not an issue in the phase-contrast image, but the fluorescence images have a glow around their outside, like reverse vignetting,^[9] that is not apparent to the eye until you turn up the contrast (or try to threshold).

With this problem and tool in mind, let me walk through the simple function FluorescentAreaMark that sits at the base of my current image processing. To get rid of the above glow I use a variant of the method we used to notice it: contrast limited adaptive histogram equalization.

def FluorescentAreaMark(img):
    #1: contrast limited adaptive histogram equalization to get rid of glow
    clahe = cv2.createCLAHE(clipLimit=50.0, tileGridSize=(20,20))

    clahe_img = clahe.apply(img)

mCherry channel after processing with adaptive histogram equalization. On the left is the image with blue corresponding to low intensity and red to high intensity, and on the right is the histogram of intensities with the cumulative distribution in black.

The algorithm works by dividing the image into 400 tiles (20 tiles by 20 tiles), and normalizing the contrast in each one. Each pixel is then given a value by bilinear interpolation between the 4 tiles nearest it. This gives a sharp distinction of the cells from their background, and we can then detect them with a simple threshold, followed by some noise canceling^[10]:

    #2: threshold and clean out salt-pepper noise
    ret, thresh_img = cv2.threshold(clahe_img,127,255, cv2.THRESH_BINARY)

    kernel = np.ones((3,3),np.uint8)

    clean_img = cv2.morphologyEx(thresh_img, cv2.MORPH_OPEN, kernel, iterations = 2)
    
    return clean_img

The threshold in line 11 simply converts the 8bit image to a binary one, by setting all pixels with intensity lower than 127 in the original image to 0 in the final image and those with intensity of 127 and higher to 1. The noise filtering step in line 15 is morphological opening — a fancy name but a straightforward algorithm. Its two basic components are erosion and dilation, both require a kernel that defines a local neighbourhood. In our case the kernel in line 13 is the one-distance Moore neighbourhood. Erosion is an AND operation on the binary mask: for a focal pixel, if every pixel in its local neighbourhood is 1 in the original image than that pixel is a 1 in the final image, else it is zero. Dilation is an OR operation: if at least one pixel in the local neighbourhood of the focal pixel is 1 in the original image then the focal pixel is a 1 in the final image. Morphological opening cycles erosion followed by dilation a number of times given by iterations — two in our case. Thus, lonely pixels that are most likely due to noise are eliminated, but the overall area of the non-noise pixels is not decreased.

The returned clean_img on line 17 is then a binary two dimensional array with a ‘1’ for each fluorescent pixel and ‘0’ otherwise. If we want the total fluorescent area then we just sum this array, outside a small buffer area around the edges which tends to contain more mistakes in correcting the edge glow. If desired, we can also combine the raw phase-contrasts images with the two processed fluorescent areas for a final visualization of the distribution of fluorescence throughout the spatially structured population:

    #get the area masks
    rdFA = FluorescentAreaMark(rdImage8)
    gnFA = FluorescentAreaMark(gnImage8)
    
    ign_buf = 30 #how big of an edge do we ignore?
    rd_area = (rdFA[ign_buf:-ign_buf,ign_buf:-ign_buf] > 0)
    gn_area = (gnFA[ign_buf:-ign_buf,ign_buf:-ign_buf] > 0)
    
    #create image to save
    img_out = []
    if ret_img:
        bW = 0.85
        
        img_out = cv2.merge((
            cv2.addWeighted(grImage,bW,rdFA,1- bW,1),
            cv2.addWeighted(grImage,bW,gnFA,1 - bW,1),
            cv2.addWeighted(grImage,bW,np.zeros_like(grImage), 1 - bW, 1)))
    
    return [np.sum(rd_area),np.sum(gn_area)], img_out

Distribution of fluorescent area. Non-small cell lung cancer cells are marked with mCherry in red, and non-cancerous fibroblasts are marked GFP in green. Colours in this image are a binary mask after the processing from this post.

In the near future, I plan to use this sort of spatial distribution information to learn about operationalizing the local environment. For the next post, however, I will just look at the overall population dynamics estimated by this method and what we can learn from them. Although the basics of computer vision that I hacked together in this post are sufficient for our purposes with this project, it does violate my ambitions as a computer scientist. If there is time and interest then I will also write some future posts to refine this code to minimize categorization errors and introduce some proper cell segmentation.

Notes

1. I was also inspired by Jan Poleszczuk’s Compute Cancer blog to start posting actual code snippets to allow for easier collaboration and engagement with programmers interested in biology, or biologists looking to start programming. I recommend checking out his blog for many code snippets and tricks across various languages. For example, if you want to oppose my choice of Python in this post then check out his speed comparison of Java, C++, MATLAB, Julia or Python for cellular automatons post for ammunition. Also, let me know what you think of interspersing code within this post. If there is a lot of interest in this format then I will try to include it in more future posts.
2. I am, of course, reinventing the wheel. A lot of great packages and programs already exist for processing and analyzing microscopy data, and we even have some in-house code here at the IMO. However, I wanted to go through the motions of this as a learning exercise. This is to make sure that I understood exactly what the algorithm is doing, and to make sure that it is a good fit for the particulars of our experimental set up.
3. Having several target audiences means that the first section on microscopy is probably too basic for biologists, and the second section on computer vision is probably too basic for programmers. So feel free to jump around. As with many of my posts, this is also meant as a tutorial for future Artem, for when I return to this experiment and code after a few months and have no idea what it does, why it exists, and how to get it running. I will try to confine the more pedantic parts of this documentation to footnotes.
4. In fact, as automation continues its perpetual march forward, even the preparation step is being eliminated. Labs like the Foundry at Imperial College London are automating the whole bench biology pipeline. I wonder if this will soon allow us to program in vitro experiments from basic components and then run them in remote automated labs, much like how we currently send in silico simulations to clusters. My optimistic hope is that the communication with machines forced by automation will also encourage more and more experimental biology graduate students to learn more and more computer science. Maybe that exposure to CS will allow for a flourishing of algorithmic biology. My cynical side, however, wonders at the disconnect between programming — a skill that everyone, from any field, should pick up — and the more niche aspects of algorithmic thinking or theoretical computer science.
5. At our magnification in this post, the field of view for a single photo is about 1/5th of the diameter of the well. The microscope takes photos of three different positions — the same positions at every 2-hour time step — in each well.
6. These limitations of not segmenting single cells can be overcome with more work. I will have to return to them in future blogposts if I want to look at single-cell tracking and lineage tracing.
7. Python has packages for almost anything, and usually in many versions, so there are plenty of alternatives to OpenCV. I chose it due to its great documentation, and from chatting with Chandler Gatenbee — the IMO expert on computer vision — who also builds his code with this package. Chandler’s warning was that the hardest part of OpenCV is getting it installed.

   I use the Anaconda distribution of Python 3, and it is pretty easy to get it with the built-in package manager from binstar: conda install -c https://conda.binstar.org/menpo opencv. It is also a good idea to have a separate environment for opencv to work in, to avoid breaking dependencies. This can be done with conda create -n opencv numpy scipy scikit-learn matplotlib python=3 and then accessed later from command line (before running jyputer notebook, for example) with source activate opencv. For a slightly more detailed tutorial, see this.

   The code I discuss in this post is available on my GitHub, and the numbered lines of code correspond to the numbers in the first commit FluorescentArea.py. The images through out this post were generated from the jyputer notebook in which I hacked the code together and paint. It is straightforward to generate these images with OpenCV, but if you want the details then ask me in the comments and I will post some more code snippets.

8. The function cv2.converScaleAbs tends to clip-down the intensities, setting every pixel above 255 to 255, so if there are intensities above 255 in the 16-bit images, it is better to first renormalize (as is done in lines 30,35 by setting alpha to something other than 1) to make the brightest pixel 255. Note that this renormalization means that I cannot make physical sense of the raw intensity numbers and have to rely on contrast. I could use a consistent renormalization to allow using raw intensity numbers, but since the rest of the processing focuses on contrast, there is no reason to worry about it.
9. We are not sure what physical process causes this glow. The apparent circular symmetry makes me think it has something to do with the lens. Andriy suspects it has more to do with partial fluorescence from the well walls, and to remedy this, we are using a different type of 96 well plate in follow up experiments.