how to simulate the following scenario in mathematica

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 ran开发者_运维问答dom 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.

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.

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

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}}

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: