how to simulate the following scenario in mathematica

Question

Suppose I have n=6 distinct monomers each of which has two distinct and reactive ends. During each round of reaction, one random end unites with another random end, either elongates the monomer to a dimer or self-associates into a loop. This reaction process stops whenever no free ends are present in the system. I want to use Mma to simulate the reaction process.

I am thinking to represent the monomers as a list of strings, {'1-2', '3-4', '5-6', '7-8', '9-10', '11-12'}, then to do one round of reacion by updating the content of the list, for example either {'1-2-1', '3-4', '5-6', '7-8', '9-10', '11-12'} or {'1-2-3-4', '5-6', '7-8', '9-10', '11-12'}. But I am not able to go very far due to my programming limitation in Mma. Could anyone please help? Thanks a lot.

score 2 · Answer 1 · answered Mar 07 '11 at 23:06

It seems like it would be more natural to represent your molecules as lists rather than strings. So start with {{1,2},{3,4},{5,6}} and so on. Then open chains are just longer lists {1,2,3,4} or whatever, and have some special convention for loops such as starting with the symbol "loop". {{loop,1,2},{3,4,5,6},{7,8}} or whatever.

How detailed does your simulation actually need to be? For instance, do you actually care which monomers end up next to which, or do you only care about the statistics of the lengths of chains? In the latter case, you could greatly simplify the state of your simulation: it could, for instance, consist of a list of loop lengths (which would start empty) and a list of open chain lengths (which would start as a bunch of 1s). Then one simulation step is: pick an open chain at random; with appropriate probabilities, either turn that into a loop or combine it with another open chain.

Mathematica things you might want to look up: RandomInteger, RandomChoice; Prepend, Append, Insert, Delete, ReplacePart, Join; While (though actually some sort of "functional iteration" with, e.g., NestWhile might make for prettier code).

I do not care which monomers end up next to which; what I care is the statistics of the lengths of the loops, and of course the number of loops. Appreciate if you could please post a working example... — Qiang Li, Mar 07 '11 at 23:31
@QiangLi: A bit of politeness wouldn't go astray. Gareth was making valid points and gave good hints. And it sounds like he's actually had experience at these types of simulations. He was just trying to clarify what you actually need. — Simon, Mar 08 '11 at 00:34
@QiangLi: Also, in a simple system like this, you could probably find an analytic solution to the statistical questions that you're asking. — Simon, Mar 08 '11 at 00:42
@Simon, @Gareth, really sorry if what I commented earlier did not sound polite. I really did not intend to. I do appreciate all of your help very much. :) — Qiang Li, Mar 08 '11 at 01:17

score 2 · Accepted Answer · answered Mar 07 '11 at 23:42

Here is the set-up:

Clear[freeVertices];
freeVertices[edgeList_List] := Select[Tally[Flatten[edgeList]], #[[2]] < 2 &][[All, 1]];

ClearAll[setNew, componentsBFLS];
setNew[x_, x_] := Null;
setNew[lhs_, rhs_] := lhs := Function[Null, (#1 := #0[##]); #2, HoldFirst][lhs, rhs];

componentsBFLS[lst_List] := 
 Module[{f}, setNew @@@ Map[f, lst, {2}]; GatherBy[Tally[Flatten@lst][[All, 1]], f]];

Here is the start:

In[13]:= start = Partition[Range[12], 2]

Out[13]= {{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}}

Here are the steps:

In[51]:= steps = 
NestWhileList[Append[#, RandomSample[freeVertices[#], 2]] &, 
  start, freeVertices[#] =!= {} &]

Out[51]= {{{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}}, {{1, 
2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {5, 1}}, {{1, 
2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {5, 1}, {3, 
4}}, {{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {5, 
1}, {3, 4}, {7, 11}}, {{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 
10}, {11, 12}, {5, 1}, {3, 4}, {7, 11}, {8, 2}}, {{1, 2}, {3, 
4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {5, 1}, {3, 4}, {7, 11}, {8,
2}, {6, 10}}, {{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 
12}, {5, 1}, {3, 4}, {7, 11}, {8, 2}, {6, 10}, {9, 12}}}

Here are the connected components (cycles etc), which you can study:

In[52]:= componentsBFLS /@ steps

Out[52]= {{{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}}, {{1, 2,
5, 6}, {3, 4}, {7, 8}, {9, 10}, {11, 12}}, {{1, 2, 5, 6}, {3, 
4}, {7, 8}, {9, 10}, {11, 12}}, {{1, 2, 5, 6}, {3, 4}, {7, 8, 11, 
12}, {9, 10}}, {{1, 2, 5, 6, 7, 8, 11, 12}, {3, 4}, {9, 10}}, {{1, 
2, 5, 6, 7, 8, 9, 10, 11, 12}, {3, 4}}, {{1, 2, 5, 6, 7, 8, 9, 10, 
11, 12}, {3, 4}}}

What happens is that we treat all pairs as edges in one big graph, and add an edge randomly if both vertices have at most one connection to some other edge at the moment. At some point, the process stops. Then, we map the componentsBFLS function onto resulting graphs (representing the steps of the simulation), to get the connected components of the graphs (steps). You could use other metrics as well, of course, and write more functions which will analyze the steps for loops etc. Hope this will get you started.

this is very novel and succinct approach! BTW, what is "BFLS" short for, Breath First List Search ...? — Qiang Li, Mar 08 '11 at 01:24
@Qiang Li No, these are the first letters of names of people involved in writing this function :) .I mentioned this function also here: http://stackoverflow.com/questions/4198961/what-is-in-your-mathematica-tool-bag/4923345#4923345 — Leonid Shifrin, Mar 08 '11 at 10:01

Simon · Answer 3 · 2011-03-08T01:32:35.267

Here's a simple approach. Following the examples given in the question, I've assumed that the monomers have a prefered binding, so that only {1,2} + {3,4} -> {1,2,3,4} OR {1,2,1} + {3,4,3} is possible, but {1,2} + {3,4} -> {1,2,4,3} is not possible. The following code should be packaged up as a nice function/module once you are happy with it. If you're after statistics, then it can also probably be compiled to add some speed.

Initialize:

In[1]:= monomers=Partition[Range[12],2]
        loops={}
Out[1]= {{1,2},{3,4},{5,6},{7,8},{9,10},{11,12}}
Out[2]= {}

The loop:

In[3]:= While[monomers!={},
  choice=RandomInteger[{1,Length[monomers]},2];
  If[Equal@@choice, 
     AppendTo[loops, monomers[[choice[[1]]]]];
       monomers=Delete[monomers,choice[[1]]],
     monomers=Prepend[Delete[monomers,Transpose[{choice}]],
                      Join@@Extract[monomers,Transpose[{choice}]]]];
     Print[monomers,"\t",loops]
   ]
During evaluation of In[3]:= {{7,8,1,2},{3,4},{5,6},{9,10},{11,12}} {}
During evaluation of In[3]:= {{5,6,7,8,1,2},{3,4},{9,10},{11,12}}   {}
During evaluation of In[3]:= {{5,6,7,8,1,2},{3,4},{9,10}}   {{11,12}}
During evaluation of In[3]:= {{3,4,5,6,7,8,1,2},{9,10}} {{11,12}}
During evaluation of In[3]:= {{9,10}}   {{11,12},{3,4,5,6,7,8,1,2}}
During evaluation of In[3]:= {} {{11,12},{3,4,5,6,7,8,1,2},{9,10}}

Edit:

If the monomers can bind at both ends, you just add a option to flip on of the monomers that you join, e.g.

Do[
  choice=RandomInteger[{1,Length[monomers]},2];
  reverse=RandomChoice[{Reverse,Identity}];
  If[Equal@@choice,
    AppendTo[loops,monomers[[choice[[1]]]]];
      monomers=Delete[monomers,choice[[1]]],
    monomers=Prepend[Delete[monomers,Transpose[{choice}]],
             Join[monomers[[choice[[1]]]],reverse@monomers[[choice[[2]]]]]]];
  Print[monomers,"\t",loops],{Length[monomers]}]

{{7,8,10,9},{1,2},{3,4},{5,6},{11,12}}  {}
{{3,4,5,6},{7,8,10,9},{1,2},{11,12}}    {}
{{3,4,5,6},{7,8,10,9},{11,12}}  {{1,2}}
{{7,8,10,9},{11,12}}    {{1,2},{3,4,5,6}}
{{7,8,10,9,11,12}}  {{1,2},{3,4,5,6}}
{}  {{1,2},{3,4,5,6},{7,8,10,9,11,12}}

While the above could be performed functionally with a `NestWhile` or `NestWhileList` (or, since you know how many steps are needed, a simple `Nest`). Since you want to gather statistics on the dynamics, you probably want to `Compile` the above into a function and so the imperative style that I used should be ok. — Simon, Mar 08 '11 at 00:40
thanks a lot. Your assumption isn't what I meant, though I might have not stated this clearly. {1,2} + {3,4} -> {1,2,4,3} and {1,2} + {3,4} -> {1,2,3, 4} are both allowed and counted as distinct. Let me try to figure out the modification to the code here... — Qiang Li, Mar 08 '11 at 01:21

Sjoerd C. de Vries · Answer 4 · 2011-03-08T09:33:40.363

I see my implementation mimics Simon's closely. Reminder to self: never go to bed before posting solution...

simulatePolimerization[originalStuff_] :=
 Module[{openStuff = originalStuff, closedStuff = {}, picks},
  While[Length[openStuff] > 0,
   picks = RandomInteger[{1, Length[openStuff]}, 2];
   openStuff = If[RandomInteger[1] == 1, Reverse[#], #] & /@ openStuff;
   If[Equal @@ picks,
    (* closing *)
    AppendTo[closedStuff,Append[openStuff[[picks[[1]]]], openStuff[[picks[[1]], 1]]]];
    openStuff = Delete[openStuff, picks[[1]]],
    (* merging *)
    AppendTo[openStuff,Join[openStuff[[picks[[1]]]], openStuff[[picks[[2]]]]]];
    openStuff = Delete[openStuff, List /@ picks]
  ]
 ];
 Return[closedStuff]
]

Some results:

enter image description here

+1 for the logical subtleties in _"never go to bed before posting a solution"_ — Dr. belisarius, Mar 08 '11 at 14:59

how to simulate the following scenario in mathematica

4 Answers4

Edit: