The curious case of missing random-walking particles in the box

Question

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):

1. specie is an Association which is of the form <| particle index or identity (integer) -> {"name of particle(either nodal or lefty)", x-y position of the created particle} ..... |>

2. specieBoundaryIndex is just an association to track whether a particle in specie 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.

3. 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

1. molecules undergo brownian motion
2. we check if the particles are inside the box and outside the cell. force them to stay
3. see ligand(nodal/lefty) interaction pairs with the receptors
4. binding mechanism to determine favourable interactions out of all pairs
5. unbinding module to remove ligand from the bound receptor
6. 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

Answer 1

The problem has been finally resolved. I used some conditions to spit out the state of the system when any particle disappeared. I found that the main reason for particles to disappear was a pattern that was defined very generally in unbindingMechanics, where a more stringent pattern was needed.

Mainly the blank after Alternatives @@ successfuldisassoc had to be replaced with the particlename.

I replaced the pattern below:

membraneIndex = Replace[membraneIndex, PatternSequence[
x_Integer -> {{p1_, p2_}, 1,Alternatives @@ successfuldisassoc,
_}] :> x -> {{p1, p2}, 0, None, "type"}, 1]

with this:

membraneIndex = Replace[membraneIndex, PatternSequence[
x_Integer -> {{p1_, p2_}, 1,Alternatives @@ successfuldisassoc,
ligandname}] :> x -> {{p1, p2}, 0, None, "type"}, 1]

If the ligandname is not specified and if two different kind of particles - with the same index - are bound to the receptors then both will get replaced

Moreover, another module biasSharedReceptors was added. This module ensures that if a membrane particle has two different types of free particles in its vicinity then only a single particle gets the opportunity to interact with it. We can generate a probability to bias the system and hence the name biasSharedReceptors.

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 *)

I have made changes to the code in the question as well; these include strong checks and corrections for the identified problems. A user can now take the code in the question and hopefully run it without any issues.