Team:LMU-TUM Munich/Model
Modeling
Since our approach to tissue printing is supposed to work without a supporting scaffold, the cohesiveness of the printed medium is crucial. Biotin binds strongly to Streptavidin as well as Avidin. Thus we are not worried about the individual bonding itself but rather about the interconnectedness of the individual cells and binding .
Before diving into the details we introduce a few abbreviations to make reading and writing more easy and efficient. Since it is yet to decide weather Streptavidin/Avidin or Biotin will be attached to the cells or to the networking (NP) we refer to the binding sites of cells as CBS and to the ones of NP as NPBS.
In the worst case all CBS we dispense from the printhead would get immediately occupied with otherwise loose NP, which would leave us with lots of individual NP-coated cells. The other extreme case would be that that the cell hardly finds any to bond to. Again resulting in individual cells that are unconnected.
So we assumed that the problem is most likely a typical question of polymerization. While the bonds in our case are not covalent we decided it should be safe to assume they were, since they are very strong with a dissociation constant of 10-15 M. So our most promising discoveries after browsing pertinent literature were the Carothers equation as well as the Flory-Stockmayer theory, which is a generalization of the former.
[describe Flory-Stockmayer theory]
[describe our assumtions to make it fit]
[describe results and limitations]
After we found that the results provided are not sufficient to optimize our process we were talking to some Professors in the field of polymer chemistry, biochemistry and biological modelling. [Most Answers were discouraging because ...] Ultimately we decided to pursue the advice of Dr. Hasenauer, who suggested to simulate the problem to understand it better. First with strong assumptions, later on reducing them stepwise to get more and more accurate results.
So our first approach assumed:
- Every cell has x CBS.
- Every NP has y PBS.
- There are 100 cells and z NP.
- Every pair of CBS and PBS have the same probability to bind. Regardless of spacial distance, obstruction, [some more fancy chemical words for this?]
- At some point either all or all cell binding sites will have bond.
In a first experiment we fixed x to 30, compromising between a way higher number of binding sites, which could “consume” , and a much lower amount of cells that could actually surround it due to spacial limitations. Furthermore we fixed y to 8 to get a first impression and run experiments with logarithmically varying z.
To analyze the result we analyzed what we call a cell-graph. The cell-graph is obtained by drawing a node for each cell and connecting each node with an edge to the nodes of the respective other cells that this cell is connected to via a .
On this graph we evaluated a couple of metrics from graph theory, resulting in the plots in Figure 1
We were happy that most metrics were pointing to roughly the same sweet spot of relative polymer/cell concentration. Soon we discovered however that the solution is a rather trivial one. Since either all CBS or all PBS are bond before the simulation stops it’s quite obvious that the sweet spot is where the total combined PBS equal the combined CBS. So next we tried to incorporate spatial effects into binding probability.
Again we assume:
- Every cell has x CBS.
- Every NP has y PBS.
- There are 100 cells and z NP.
But instead of assuming equal probability for each CBS/PBS pair to bind we assume:
- A one-dimensional space, where all cells as well as NP are placed in equidistant manner.
- For now the position of individual CBS is assumed equal to the respective cell. Likewise for PBS and NP.
- The size of the space is defined by a constant c, where 1 means that the space is of size 1 and 0.01 means the space is of size 100. This constant corresponds to the concentration of active components in a medium. The lower the concentration the higher the distance between two components.
- At every timestep a constant fraction of free PBS tries to bind. (reactivity)
- When a PBS is trying to react, each free CBS gets assigned a certain probability depending on two principles:
- The more near a CBS is to the PBS, the more likely it is to bind to this one. [zero mean gaussian of the distance with a certain varaiance \sigma]
- The more cells are bound to a certain NP already, the less likely it is, that it will bind to one that it didn’t bind to yet. [limited space, big cells, small NP] [which distribution?]
- With increasing time, the variance \sigma is increasing to modell that cells and NP can move freely in the medium and reach more distant partners. [neglecting that they only move in one direction; let’s say we don’t know which one]
- After a defined numebrr of timesteps the experiment is stipped and the result is assessed the same way described under [section above]
References
Code
<!DOCTYPE html>
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from monomers import *
def init_monomers(n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4):
Cell.n_cells=0
cells=[]
for i in range(n_cells):
cells.append(Cell(n_cell_binding))
.n_proteins=0
=[]
for i in range(n_proteins):
.append((n_protein_binding))
return cells,
def polymerize(cells, ):
free_cell_bindings = []
for c in cells:
free_cell_bindings += c.free_bindings
free_protein_bindings = []
for p in :
free_protein_bindings += p.free_bindings
while len(free_protein_bindings)>0 and len(free_cell_bindings)>0:
p_binding_idx = np.random.randint(len(free_protein_bindings))
p_binding = free_protein_bindings.pop(p_binding_idx)
c_binding_idx = np.random.randint(len(free_cell_bindings))
c_binding = free_cell_bindings.pop(c_binding_idx)
#print c_binding.parent.id, p_binding.parent.id
p_binding.bind(c_binding)
return len(free_protein_bindings) , len(free_cell_bindings)
def polymerize_spacial(cells, , n_protein_binding, concentration=0.01, reactivity=0.2, mobility=0.05, max_time=10.):
time = 0.
position_variance = 1.
distance_binding_probability = lambda x: np.exp(-x ** 2 / position_variance ** 2) / (
2 * np.pi * position_variance ** 2)
free_protein_bindings = []
for p in :
free_protein_bindings += p.free_bindings
while (time < max_time):
for i in range(int(np.ceil(reactivity * len(free_protein_bindings)))):
# choose a random PBS
p_binding_idx = np.random.randint(len(free_protein_bindings))
p_binding = free_protein_bindings[p_binding_idx]
# do tests on parent
p = p_binding.parent
# position: the liwer to concentration, the higher the distance
p_pos = float(p.id) / (len() * concentration)
# (float) idx of nearest cell
c_id_mean = p_pos * len(cells) * concentration
# max distance of reachable cells from nearest cell
max_range = position_variance * 3.2
c_id_range = max_range * len(cells) * concentration
c_id_min = max(int(c_id_mean - c_id_range), 0)
c_id_max = min(int(c_id_mean + c_id_range), len(cells))
candidate_cells = cells[c_id_min:c_id_max]
# each has a probability according to the distance
if len(candidate_cells) == 0:
continue
min_distance = max(p_pos - position_variance * 3.2, 0) - p_pos
max_distance = min(-p_pos + position_variance * 3.2,
len() - 1 / len() / concentration) + p_pos
probabilities = np.arange(min_distance, max_distance,
(max_distance - min_distance) / len(candidate_cells))[:len(candidate_cells)]
probabilities = distance_binding_probability(probabilities)
probabilities /= sum(probabilities)
# pick one
cell = np.random.choice(candidate_cells, p=probabilities)
try:
cell.free_bindings.next()
except StopIteration:
# no free bindings left in this cell
continue
# let react with probability respective to free bindings
protein_partners = list(set([b.partner.parent for b in p.bindings if b.partner is not None]))
if len(protein_partners) == 1:
# if a is already attatched to a cell it will most likely not bind to other cells
if cell not in protein_partners and np.random.rand()<0.9:
continue
if len(protein_partners)>1 and cell not in protein_partners:
# if a was binding to two cells already, it will not bind to any other
continue
# if np.random.rand() < float(4 - len(list(p.free_bindings))) / 3.:
# # can only bind to cells it's already bound to
# if cell not in [b.partner.parent for b in p.bindings if b.partner is not None]:
# continue
cell.free_bindings.next().bind(p_binding)
free_protein_bindings.pop(p_binding_idx)
position_variance += mobility
time += 0.01
free_cell_bindings = []
for c in cells:
free_cell_bindings += c.free_bindings
return len(free_protein_bindings), len(free_cell_bindings)
draw = True
_,_,_=random_graph(cell_graph, 30, 120, n_cell_binding=400, n_protein_binding=4, concentration=0.02, reactivity=0.1, mobility=0.05, max_time=5.0)
Analysis¶
draw = False
from graph_tool.all import *
from graph_tool import clustering,stats
def cell_graph(cells,,draw=False):
'''vertices are cells, cells connected via a are connected in the graph with an edge'''
g = Graph(directed=False)
vertices = list(g.add_vertex(len(cells)))
for c_id in range(len(cells)):
for c_b in cells[c_id].bindings:
try:
for p_b in c_b.partner.parent.bindings:
try:
c2_id = p_b.partner.parent.id
if not c_id==c2_id:
g.add_edge(vertices[c_id],vertices[c2_id])
except AttributeError:
pass
except AttributeError:
pass
if draw:
graph_draw(g, vertex_text=g.vertex_index, edge_pen_width=2 , vertex_font_size=30, output="/tmp/graph.png",
output_size=(2000, 2000))
return g
#graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
# output_size=(2000, 2000))
def random_graph(graph_method, n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4,concentration=0.02,reactivity=0.01,mobility=0.1,max_time=2.):
# cells,=init_monomers(n_cells, n_proteins, n_cell_binding, n_protein_binding)
# polymerize(cells,)
cells,=init_monomers(n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4)
polymerize_spacial(cells, , n_protein_binding=n_protein_binding, concentration=concentration, reactivity=reactivity, mobility=mobility, max_time=max_time)
return graph_method(cells,,draw),cells,
results=[]
draw=False
for n_proteins in np.arange(-1.5,1.5,0.1):
n_proteins = int(len(cells)*10**n_proteins)
g,cells, = random_graph(cell_graph, 30, n_proteins, n_cell_binding=400, n_protein_binding=4, concentration=0.02, reactivity=0.1, mobility=0.05, max_time=5.0)
global_clustering = clustering.global_clustering(g)[0]
tmp=graph_tool.topology.max_cardinality_matching(g)
n_max_matching = tmp.count(tmp)
max_shortest_path = graph_tool.topology.pseudo_diameter(g)[0]
tmp=graph_tool.topology.label_largest_component(g)
biggest_cluster = tmp.get_array().shape[0]
tmp = graph_tool.topology.label_components(g)[1]
tmp = np.sort(tmp)[::-1]
n_clusters = len(tmp)
biggest_cluster = tmp[0]
mean_cluster = np.mean(tmp)
n_free_protein_bindings = 0
n_free_proteins = 0
for p in :
p_free = len(list(p.free_bindings))
n_free_protein_bindings += p_free
if p_free == len(p.bindings):
n_free_proteins += 1
n_free_cell_bindings = 0
n_free_cells = 0
for c in cells:
c_free = len(list(c.free_bindings))
n_free_cell_bindings += c_free
if c_free == len(c.bindings):
n_free_cells += 1
print n_proteins, n_free_cells, n_free_cell_bindings, n_free_proteins, n_free_protein_bindings, n_clusters, mean_cluster, biggest_cluster, global_clustering, max_shortest_path#, n_max_matching, biggest_cluster
results.append([n_proteins, n_free_cells, n_free_cell_bindings, n_free_proteins, n_free_protein_bindings, n_clusters, mean_cluster, biggest_cluster, global_clustering, max_shortest_path])
tmp_res=results
results=np.asarray(results).T
n_proteins=results[0]
results=results[1:]
titles=["n_free_cells", "n_free_cell_bindings", "n_free_proteins", "n_free_protein_bindings", "n_clusters", "mean_cluster", "biggest_cluster", "global_clustering", "max_shortest_path"]
f, axes = plt.subplots(len(titles)/2+1,2, figsize=(16,8*(len(titles)/2+1)))
for i,result in enumerate(results):
ax = axes[i/2,i%2]
ax.set_title(titles[i])
ax.plot(n_proteins,result.T,"-x")
ax.set_xscale("log")
clustering.global_clustering(g)
stats.distance_histogram(g)
tmp=graph_tool.topology.max_cardinality_matching(g)
tmp.count(tmp)
graph_tool.topology.max_independent_vertex_set(g)
graph_tool.topology.pseudo_diameter(g)
tmp=graph_tool.topology.label_largest_component(g)
tmp.get_array().shape[0]
Cells and bindings vertices¶
g = Graph(directed=False)
v_size = g.new_vertex_property("int",val=20)
v_group = g.new_vertex_property("int32_t")
#v_color = vprop_vint = g.new_vertex_property("<int>", vals=[])
#cell_vertices = list(g.add_vertex(n_cells))
for cell in cells:
cell.vertex = g.add_vertex()
v_size[cell.vertex] = 100
for p in :
last_vertex = None
for b in p.bindings:
b.vertex=g.add_vertex()
v_size[b.vertex] = 10
v_group[cell.vertex] = p.id
#intraprotein connections
if last_vertex is not None:
g.add_edge(b.vertex,last_vertex)
last_vertex=b.vertex
#cell connections
if b.partner is not None:
cell_vertex = b.partner.parent.vertex
g.add_edge(b.vertex,cell_vertex)
pos=sfdp_layout(g,groups=v_group)
#graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
# output_size=(2000, 2000))
graph_draw(g, pos=pos,edge_pen_width=2, vertex_text=g.vertex_index, vertex_size=v_size, vertex_font_size=0, output_size=(2000, 2000))
#graphviz_draw(g, elen=100, vsize=v_size, size=(2000, 2000))
vertices and cell vertices¶
def pc_graph(cells, , draw=False):
g = Graph(directed=False)
v_size = g.new_vertex_property("int")
#v_group = g.new_vertex_property("i'nt32_t")
#v_color = vprop_vint = g.new_vertex_property("<int>", vals=[])
for cell in cells:
cell.vertex = g.add_vertex()
v_size[cell.vertex] = 120
for p in :
p.vertex=g.add_vertex()
v_size[p.vertex] = 40
for b in p.bindings:
#cell connections
if b.partner is not None:
cell_vertex = b.partner.parent.vertex
if g.edge(p.vertex,cell_vertex) is None: #only one connection per pair
g.add_edge(p.vertex,cell_vertex)
pos=fruchterman_reingold_layout(g)
#graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
# output_size=(2000, 2000))
if draw:
graph_draw(g, edge_pen_width=2, vertex_text=g.vertex_index, vertex_size=v_size, vertex_font_size=15, output_size=(2000, 2000))
return g
g.vertex_index[vertices[1]]
cells[23].bindings[5].partner.parent.bindings[0].partner.parent
