Note: I added the entire code I wrote below to generate the initial configuration and have explained step by step what each part of the code does
Kindly bear with the code and the explanation. Since the code may be dense, I will make my utmost effort here to summarize the problem that has stumped me for the past three days even after carefully reviewing my code multiple times.
I am trying to simulate a reaction-diffusion system in which there are two particles that can roam freely (call them ligands) shown in blue and red (in the picture below). Particles on the membrane are called receptors (purple). The ligands are present within a virtual biological cell (shown by a circle) and they diffuse outward from the cell. The ligands once outside the cell can bind to the membrane receptors. In fact the red particle (ligand known as lefty in biology) competes with the blue particle (known as nodal) for the same receptor (purple).
I create the initial configuration for blue and red using initialParticleConfig
function and for the purple receptors I use membraneReceptorConfig
. All code will be present at the end with comments and bold title on top declaring what the block does.
The rules of the game: When red(lefty) binds to a particular purple dot then that purple dot becomes unavailable temporarily for binding with any other particle. It is similar for the case of blue dot as well i.e. when it binds to purple dot then the purple dot temporarily becomes unavailable. I am using the word temporarily because the red and blue dots can get released from the already bound purple dots making them available again
DATA STRUCTURES
the data structures i have used are pretty straightforward. initialParticleConfig
yields 3 things at the end when provided with the appropriate conditions (number of particles and where to create them):
specie
is anAssociation
which is of the form<| particle index or identity (integer) -> {"name of particle(either nodal or lefty)", x-y position of the created particle} ..... |>
specieBoundaryIndex
is just an association to track whether a particle inspecie
is inside the cell or outside the cell. Therefore it is of the form <|particle index (an integer) -> 0 ||1 ..... |>. The value 0 is when a particle is inside the cell and 1 when it is outside.speciePts
is the list of the x-y coordinates of the positions of the point.
we first use the same code block initialParticleConfig
to generate both red and blue dots
membraneReceptorConfig
also takes in geometrical parameters and a number (specifying how many receptors to create) to place RandomPoints
on the circle. The function in general returns two things: list of all the receptor positions created on the cell (circle) and a second structure which is a list of rules where the rule corresponds to:
membrane receptor index -> {x-y coordinates of receptor, 0 || 1, None, "type"}
. The 0
here represents that the receptor is freely available to bind with red or blue dot. 1
represents that the receptor is already in use. At the start of the simulation all receptors are free. None
can be replaced with the particle-index of the blue or red dot to which the receptor is bound to. "type" can be either "nodal" or "lefty" (the names of the particles).
membranePts
delineate the position of the receptors on the membrane. Receptors remain fixed. they do not move
The results from the initialParticleConfig
and membraneReceptorConfig
are stored in global variables, namely nodallist
, nodalBoundaryindex
, nodalpts
for nodal and likewise for lefty. membranepts
, membraneindex
Strategy for the code
- molecules undergo brownian motion
- we check if the particles are inside the box and outside the cell. force them to stay
- see ligand(nodal/lefty) interaction pairs with the receptors
- binding mechanism to determine favourable interactions out of all pairs
- unbinding module to remove ligand from the bound receptor
- repeat the process
Now THE PROBLEM
With simple diffusion and checks to keep the particles in the box, nothing strange happens. However, the moment when I try to simulate my particle with the unbinding module on i.e. not commented, I tend to lose particles over time. I have seen my code more than 10 times but cannot seem to identify problem.
Please see the image below. The initial number of nodal I started from is 54 but at some time when i stop the simulation (Manually closing the Animate cell) and counted the total nodal in the system, it came out to be 53. I have tried many times and i see a decrease in either nodal or lefty or both. Please note there is no destruction term for the particle in my system. Therefore I expect the total number to remain conserved.
Also if i uncomment the unbinding part and binding part the ligands behave as expected they stay inside the box and remain outside the cell (once they diffuse out of it). And the binding part seems to work ok too. So I do not know where in the unbinding part or elsewhere lies the problem
Code that will be used to generate intial Configuration
(* setting up initial configuration of ligands *)
initialParticleConfig[name_String, {x_?NumberQ, y_?NumberQ},
number_?IntegerQ] :=
Module[{specie, specieBoundaryIndex, speciePts},
specie = RandomReal[{x, y}, {number, 2}]; (*this creates a number of
particles for specified domains *)
specie = Select[specie, # \[Element] \[ScriptCapitalR] &]; (*
select those particles that are present within the cell *)
specie = Association@MapIndexed[First@#2 -> #1 &, Thread[{name, specie}]];
(*lets thread the particles with a rule where the particles are
associated with an index, their name and their coordinates *)
specieBoundaryIndex = (#[[2]] \[Element] \[ScriptCapitalR]) & /@
Values[specie] /. {True -> 0, False -> 1}; (*a mapping over particle
position to check whether the particles are outside the defined boundaries
(cells) or within it. a particle inside would yield True (replace True with
1) and external to the region will yield False (replace with 0) *)
specieBoundaryIndex = Association@Thread[Range[Length@specie] ->
specieBoundaryIndex]; (*putting an index to the previous result for unique
identification of particle *)
speciePts = Values@specie[[All, 2]]; (*this will only extract the particle
positions *)
Return@{specie, specieBoundaryIndex, speciePts}; (*now i can return three
variables. specie (particle number \[Rule] {{"Nodal || Lefty"},{x,y}} ...),
specieBoundaryIndex (particle number \[Rule] 1 || 0 ...) and
speciePts ({x,y}....) *)
];
(* setting up configuration of membrane proteins *)
membraneReceptorConfig[{x_?NumberQ, y_?NumberQ}, radius_?NumberQ,
number_?IntegerQ] := Module[{membranePts, membraneIndex},
membranePts = Table[RandomPoint[Circle[{x, y}, radius]], number]; (*this
creates a random number of particle on a circle, basically representing
membrane proteins *)
membraneIndex = Thread[Range[number] ->
Thread[{membranePts, 0, None, "type"}]];
(* putting an identification number for the previous result i.e. receptor
number \[Rule] {coord, 0 or 1, BoundTo/None, Type-Ligand} *)
Return@{membranePts, membraneIndex} ;
(* returns membranePts ({x,y}....) and membraneIndex
(membrane protein number \[Rule] {x,y}...) *);
];
Region Constraint
\[ScriptCapitalR] = ImplicitRegion[x^2 + y^2 <= 0.25^2, {x,
y}]; (* defining the region containing the molecules in a single cell *)
Subscript[\[ScriptCapitalR], 2] =
ImplicitRegion[x^2 + y^2 <= 0.75^2, {x, y}];
Subscript[\[ScriptCapitalR], 3] =
BoundaryMeshRegion[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}}, Line[{1, 2, 3,
4,1}]];
Generate Initial Configuration
{nodallist, nodalBoundaryindex, nodalpts} =
initialParticleConfig["nodal", {-0.25, 0.25}, 70];
{leftylist, leftyBoundaryindex, leftypts} =
initialParticleConfig["lefty", {-0.25, 0.25}, 50];
{membranepts, membraneindex} =
membraneReceptorConfig[{0, 0}, 0.75, 30];
Graphics[{{Gray, Dashed,
Line[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}]}, {Black,
Opacity[0.2], Circle[{0, 0}, 0.75]}, {PointSize[Medium],
Opacity[0.6], Darker@Purple, Point@membranepts}, {Darker@Blue,
PointSize[Medium], Opacity[0.6], Point@nodalpts}, {Red,
Opacity[0.6], PointSize[Medium], Point@leftypts}},
ImageSize -> Large]
THE CODE
once the global variables or the initial configurations are created together with regions (Circle and Rectangular box to which everything is confined), a series of operations are performed to get the final state of the system.
The following line of code is used to run the simulation:
code to start simulation
tracknodal = {}; tracklefty = {};
Monitor[For[i = 1, i < 20000, i++,
g = BrownianSimulation[nodallist, nodalBoundaryindex, leftylist,
leftyBoundaryindex, membraneindex, membranepts];
], g]
I will mention the main first and all the individual operations will be listed later. Note the local variables are for most part storing values of the global variables:
BrownianSimulation[nodalp_, nodalBoundaryIndex_, leftyp_,
leftyBoundaryIndex_, membraneind_, membranep_] :=
Module[{nodalcopy, nodal, lefty, leftycopy, nodalthreadlist,
leftythreadlist, interactingNodalReceptor, interactingLeftyReceptor,
leftyind, nodalind, membraneIndex, preNodcount, newNodcount,
newpoints, newindices, newptslist, boundnodal, boundlefty,
difference},
nodalind = nodalBoundaryIndex; (*initial boundaryconfiguration of nodal *)
leftyind = leftyBoundaryIndex; (*initial configuration of lefty*)
nodal = nodalp; (* save initial configuration of nodal*)
lefty = leftyp; (* likewise *)
membraneIndex = membraneind; (*lets store configuration of membrane
receptors*)
nodalcopy = nodal; (* a copy that we can use to restore the position if the
step the particle takes is outside the box or if it wants to come to the
inside of the circle *)
leftycopy = lefty; (*likewise*)
(* ------ brownian step: the function takes in free particles and adds
randomWalk to them ------ *)
{nodal, lefty} = Map[If[Length@# != 0, BrownianCheck[#], #] &,
{nodal, lefty}];
(* ------ ligand spatial checks: this makes sure that the particles outside
the cell stay outside and also not move out of the square ------ *)
{nodalthreadlist, leftythreadlist} =
MapThread[ligandCheck, {{nodal, lefty}, {nodalcopy, leftycopy}, {nodalind,
leftyind}}];
nodalind = Association@Map[#1[[1]] -> #[[2]] &, nodalthreadlist];
leftyind = Association@Map[#1[[1]] -> #[[2]] &, leftythreadlist];
(* ------ check for possible ligand-receptor interactions: what we do here is
to find all the free receptors and the ligand (nodal or lefty) and see what
possible interactions exist. We use the NearestFunction for this. Note: an
individual receptor can only interact with one ligand ------ *)
{interactingNodalReceptor, interactingLeftyReceptor} =
Map[(nearestFunction[#, membraneIndex, membranep] /.
nearestFunction[__] :>{}) &, {nodalthreadlist,leftythreadlist}];
(* biasing a common receptor interacting with nodal and lefty at the same
time *)
{interactingNodalReceptor, interactingLeftyReceptor} =
biasSharedReceptors[interactingNodalReceptor,interactingLeftyReceptor] /.
biasSharedReceptors[pattern__] :> {pattern};
(* ------ binding mechanics: here we take the previous interactions
calculated directly above to see which interactions are favourable i.e.
we find the ligand-receptor pair which bind successfully. For this we
generate probabilities for each pair. If a successful interaction is found
then the membraneIndex (receptor is connected with the particular
particle, itd availability is altered to 1 which means not available for more
interactions), and nodalthreadlist, and nodalind are updated (particle index
are removed.the data stays preserved because we have connected it with the
membraneIndex)
Also Note: if no interactions were observed i.e.
interactingNodalReceptor or interactingLeftyReceptor is empty then the rule
below will preserve the parameters being passed ------ *)
{nodalthreadlist, nodalind, membraneIndex} = (interAct[##] & @@
{nodalthreadlist, nodalind, membraneIndex,
interactingNodalReceptor, "nodal"}) /.interAct[patt__, y_, z_] :> {patt};
{leftythreadlist, leftyind,membraneIndex} = (interAct[##] & @@
{leftythreadlist, leftyind, membraneIndex,
interactingLeftyReceptor,"lefty"}) /.interAct[patt__, y_, z_] :> {patt};
(* ------ unbinding of particles: in this we find all the bound particles in
the system and generate for each of them a probability to unbind. If the
unbinding process is a success add the ligand/particle back from the
membraneIndex (the receptor to which it was bound) to nodalthreadlist or
leftythreadlist, and likewise for nodalind. For membraneIndex update the
{availability from 1 to 0, particle index to None,"nodal"||"lefty" to "type"}
----- *)
{membraneIndex, nodalthreadlist,nodalind} =
(unbindingMechanics[membraneIndex, nodalthreadlist, nodalind,
"nodal"]) /. unbindingMechanics[patt__, y_] :> {patt};
{membraneIndex, leftythreadlist,leftyind} =
(unbindingMechanics[membraneIndex, leftythreadlist,leftyind,
"lefty"]) /. unbindingMechanics[patt__, y_] :> {patt};
(* system stats *)
(* finding bound nodal and lefty from membraneIndex *)
boundnodal = Cases[membraneIndex,
patt : PatternSequence[_ -> {{_, _}, _, _,"nodal"}] :> patt[[2, 1]], 2];
boundlefty = Cases[membraneIndex,
patt : PatternSequence[_ -> {{_, _}, _, _,"lefty"}] :> patt[[2, 1]], 2];
AppendTo[tracknodal, Length@nodalthreadlist];
AppendTo[tracklefty, Length@leftythreadlist];
(* lets change the nodallist and leftylist back to Association so that the
inputs to the BrownianSimulation remain the same *)
nodallist = Association@Map[#[[1]] -> {"nodal", #[[3]]} &, nodalthreadlist];
leftylist = Association@Map[#[[1]] -> {"lefty", #[[3]]} &, leftythreadlist];
nodalBoundaryindex = nodalind;
leftyBoundaryindex = leftyind;
membraneindex = membraneIndex;
membranepts = membraneIndex[[All, 2, 1]];
(* drawing the system state *)
Graphics[{{Gray, Dashed,
Line[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}]}, {Black,
Opacity[0.2], Circle[{0, 0}, 0.75]}, {PointSize[Medium],
Opacity[0.6], Darker@Purple,
Point@membraneIndex[[All, 2, 1]]}, {Darker@Blue,
PointSize[Medium], Opacity[0.6],
Point@nodalthreadlist[[All, 3]]}, {Red, Opacity[0.6],
PointSize[Medium],
Point@leftythreadlist[[All, 3]]}, {Darker@Green,
PointSize[Large], Opacity[0.4], Point@boundnodal}, {Darker@Red,
PointSize[Large], Opacity[0.4], Point@boundlefty}},
ImageSize -> Large]
]
first operation: Brownian walk to make all particles move
BrownianWalk[number_, elem_] :=
Module[{normal, oldposition = elem[[2]], newposition, step,
tag = elem[[1]]},
step = RandomVariate[NormalDistribution[0, 0.01], 2];
newposition = oldposition + step;
number -> {tag, newposition}
]
BrownianCheck[assoc_] := Module[{temp = assoc},
temp = Association[KeyValueMap[BrownianWalk, temp]]
] /; Length@assoc > 0
second operation: To keep particls inside the box and outside the cell
ligandCheck[particles_, speciecopy_, particleBoundaryIndex_] :=
Module[{specie = particles, speciekeys, coordsspecie,indexspecie,threadList},
speciekeys = Keys@specie;
coordsspecie = Values@specie[[All, 2]];
indexspecie = Values@particleBoundaryIndex;
threadList = Thread[{speciekeys, indexspecie, coordsspecie}];
threadList = threadList /. {ind_?IntegerQ, x_?IntegerQ, y_List} /;x == 0 &&
!y \[Element] Subscript[\[ScriptCapitalR], 2] :> {ind, 1, y};
threadList =
threadList /. {ind_?IntegerQ, x_?IntegerQ, y_List} /;
x == 1 && y \[Element] Subscript[\[ScriptCapitalR], 2] :>
{ind,1, speciecopy[ind][[2]]};
threadList =
threadList /. {ind_?IntegerQ, x_?IntegerQ, y_List} /;
x == 1 && !y \[Element] Subscript[\[ScriptCapitalR], 3] :>
{ind, 1,speciecopy[ind][[2]]};
threadList
]
third operation: nearest function to check interactions between receptors and pairs
nearestFunction[ligandList_, membraneIndex_, membranepts_] :=
Module[{freeReceptors, receptormap, interactions, ligandlist},
(* ligands *)
ligandlist = Cases[ligandList, {ind_, 1, pos_} :> {ind, pos}];
freeReceptors = Cases[membraneIndex,
PatternSequence[x_ -> {y_List, 0, None, "type"}] :> y -> x] ; (*finding
all the free receptors that are available for binding and making a
receptor map for nearest function *)
receptormap = Nearest[freeReceptors]; (* creates a nearest function for
receptors *)
interactions =
Reap[Map[
Function[{particle},
Sow[particle[[1]], receptormap[particle[[2]], {1, 0.01}]]],
ligandlist], _, List][[2]] //
DeleteDuplicatesBy[Replace[#,
{x_Integer, {y_Integer, z___Integer}} :> {x,y},2],Last] &;
Reverse[interactions, 2]
]/;Count[membraneIndex, 0, 3] > 0
before fourth operation: biasSharedReceptors: if a receptor is shared by two particles then favour a single interaction
biasSharedReceptors[nod_, left_] :=
Module[{nodal = nod, lefty = left, union, probability},
union = Cases[GatherBy[Join[nodal, lefty], Last], {_, x_} ..];
Map[
Function[{part},
probability = RandomReal[];
If[probability > 0.5,
lefty =
Delete[#,
Flatten@Position[#, Flatten@Intersection[part, #]]] &@lefty,
nodal =
Delete[#,
Flatten@Position[#, Flatten@Intersection[part, #]]] &@nodal]
], union];
Return@{nodal, lefty};
] /; (nod =!= {} && left =!= {}) (* If a receptor is interacting with
both nodal and lefty at the same time, bias the system to choose one *)
fourth operation: to determine favourable interactions between receptor-particle pairs
bindingEvent[membraneIndex_, {ligandpos_, receptorpos_}, ligandnom_] :=
Module[{coord},
coord = membraneIndex[[receptorpos]][[2, 1]];
ReplacePart[membraneIndex,
receptorpos -> (receptorpos -> {coord, 1, ligandpos, ligandnom})]
]
interAct[specie_, specieind_, membraneind_,
interactingspecies : {{_, _} ..}, ligandname_] :=
Module[{interactionProb, successpos, successInt, membraneIndUpdate ,
successligand, specielist = specie, index = specieind},
interactionProb = Table[RandomReal[], Length@interactingspecies];
successpos = Position[interactionProb, _?(# > 0.5 &)];
successInt = interactingspecies[[Flatten@successpos]];
membraneIndUpdate =
Fold[bindingEvent[#1, #2, ligandname] &, {membraneind, ## & @@successInt}];
successligand = # & @@@ successInt;
specielist = Fold[Function[{list, index},
list /. {y_, _, {_, _}} /; y == index :> Sequence[]],
{specielist, ## & @@ successligand}];
KeyDropFrom[index, successligand];
{specielist, index, membraneIndUpdate}
]/; Count[membraneind, 0, 3] > 0
last operation: to unbind any receptor from particle
unbindingMechanics[membraneind_, ligandlist_, ligandind_,
ligandname_?StringQ] :=
Module[{membraneIndex = membraneind, ligandList = ligandlist,
ligandInd = ligandind, boundSpecies, probability, position,
successfuldisassoc},
boundSpecies = Cases[membraneIndex,
PatternSequence[_Integer -> {{p1_, p2_}, 1, y_,ligandname}] :> {y, {p1,p2}}];
probability = Table[RandomReal[], Length@boundSpecies];
position = Position[probability, _?(# > 0.97 &)];
successfuldisassoc = boundSpecies[[All, 1]][[Flatten@position]];
membraneIndex = Replace[membraneIndex, PatternSequence[
x_Integer -> {{p1_, p2_}, 1,Alternatives @@ successfuldisassoc,
ligandname}] :> x -> {{p1, p2}, 0, None, "type"}, 1];
ligandInd = KeySort@AssociateTo[ligandInd, Thread[successfuldisassoc -> 1]];
ligandList = SortBy[Fold[Insert[#1, {#2[[1]], 0, #2[[2]]}, -1] &,
ligandList, boundSpecies[[Flatten@position]]] , First];
{membraneIndex, ligandList, ligandInd}
]/; Count[membraneind, ligandname, 3] > 0
tracknodal
and lefty seem not to be defined anywhere. $\endgroup$tracknodal ={}
andtracklefty={}
in the beginning of the program. thanks ! $\endgroup$nearestFunction
. Do you think a free particle replacing an already bound receptor might be destroying the bound particle? However, it is quite strange that it happens because I always make my nearestFunction from the free receptors and not the unavailable ones i.e.Nearest[freereceptors]
$\endgroup$