# YASARA MACRO # TOPIC: 4. Molecular modeling # TITLE: Packing library # REQUIRES: Dynamics # AUTHOR: Elmar Krieger & Kornel Ozvoldik # LICENSE: GPL # DESCRIPTION: This macro provides a function library for assembling giant structures with millions to billions of atoms, to be 'included' by other macros # Include library functions include md_library # BUILD A TRANSPARENT BALL&STICK MODEL # ==================================== # Build a mesh model consisting of spheres and cylinders in place of atoms and bonds. # Select atoms with 'selection' and model transparency by the blending factor 'alpha'. # Sphere radii will be equal to atom radii times 'atomscale', # cylinder radii will be equal to 'bondradius'. def BuildBallSticks selection,alpha,atomscale,bondradius # Convert selection to atom list atomlist()=ListAtom (selection) # Loop over the atoms for atom in atomlist # Create atom sphere radius=RadiusAtom (atom) color=ColorAtom (atom) pos()=PosAtom (atom),CoordSys=Global obj=ShowSphere Radius=(radius*atomscale),Color=(color),Alpha=(alpha),Level=2 # Shrink radius to ensure that bond cylinders touch the atom spheres radius=sqr (radius*atomscale)-sqr bondradius if radius<0 radius=0 else radius=sqrt (sqr (radius*atomscale)-sqr bondradius) PosObj (obj),(pos) # Loop over the bonds linkedlist()=ListAtom all with bond to (atom) for linked in linkedlist if linked1 # There is more than one position, orient filament along first path endpos()=poslist4,poslist5,poslist6 ori()=OriVec (endpos-pos) # Each two positions define one path paths=positions-1 # Set starting position and orientation PosOriObj (startobj),(pos),(ori) # Replace all duplicate molecules by instances MonomerizeObj (startobj) OligomerizeObj (startobj),Instance=Yes instances=InstanceObj (startobj) for i=1 to paths if positions>1 # Build filament along paths defined by positions in 'poslist' for j=1 to 3 # Get end position of current path endpos(j)=0.+poslist(i*3+j) # Get direction vector to end position of current path newdir()=endpos-pos # Calculate how many fragments can be placed along current path dirlen=norm newdir fragments=0+round (dirlen/steplen) # Normalize direction vector newdir()=newdir()/dirlen else # Calculate how many fragments we need fragments=1+round ((0.+length)/lenperfrag) newdir1,newdir2,newdir3=VecOri (ori) for j=1 to fragments dir()=newdir # Move ahead pos()=pos+dir*steplen if positions>1 # Get direction vector to end position of current path if j==fragments k=0 if isteplen/2) break if i==paths or k==paths # For the last fragment of the last path keep the same direction newdir()=dir else # Apply normal rotation around X-axis ori1=ori1+alpha # Apply random rotation ori()=ori+rotintstart+rnd rotintlen # New direction vector with current orientation 'ori' newdir1,newdir2,newdir3=VecOri (ori) # Rotation axis is given by the cross product of the old and new dir vector axis()=newdir(2)*dir(3)-newdir(3)*dir(2),newdir(3)*dir(1)-newdir(1)*dir(3),newdir(1)*dir(2)-newdir(2)*dir(1) # Rotation angle is given by the inner product of the old and new dir vector angle=-acos (sum (newdir*dir)) # Get position and orientation of last set of instances totinstances=InstanceObj (startobj) posorilist()=PosOriObj (startobj) Ins (1+totinstances-instances)-(totinstances) # Create new set of instances InstanceObj (startobj),(instances) for k=1 to instances ins=totinstances+k # Set position and orientation for each instance for l=1 to 6 posori(l)=posorilist((k-1)*6+l) PosOriObj (startobj) Ins (ins),(posori) # Move and rotate instance MoveObj (startobj) Ins (ins),(dir*steplen) RotAxisObj (startobj) Ins (ins),(pos),(dir),(alpha) if norm axis RotAxisObj (startobj) Ins (ins),(pos),(axis),(angle) NameObj (startobj),(filament) HideMessage return (startobj) # INJECT A GENOME INTO A GROWING MEMBRANE BUD # =========================================== # To create a realistically folded genome, it is shifted into a growing membrane bud # during a coarse grained MD simulation. 'radius' is the final radius of the bud. # The genome must have been created with BuildGenome and no other objects may be active. # If 'turned', the genome is turned around (in case NiceOri chooses the wrong direction) def InjectGenome radius,turned global Pi # Orient genome along the X-axis NiceOriAll Pos 0,0,50 if turned Rotate X=180 Rotate Y=180 # Get the width and position along X, place the genome start at the center genlenx,_,_,genposx=GroupBox all,CoordSys=Global MoveAll (-genposx+genlenx/2) # Create a cell large enough for the bud cubelen=2.*radius+10 Cell (cubelen) Boundary Wall MoveObj SimCell,X=(-cubelen/2) plane=ShowBox Width=(cubelen*2),Height=10,Depth=(cubelen*3),LeftCol=Gray MoveObj (plane),X=(-cubelen/8),Y=(-cubelen/2-5) LightSource CellShadow=Yes # The bud eccentricity, it is moved a bit off center so that the injected genome does to bump straight into the other side budecc=4 # Determine tesselation level of the sphere that represents the membrane bud. # The larger 'radius', the more vertices (and ultimatly atoms) do we need to make sure # that the genome cannot escape through holes when the sphere reaches its final size. sphere=ShowSphere (radius),Red,Level=(0+log (sqr radius)*0.6) # Turn the polygon mesh into atoms that can serve as a boundary during the MD simulation bud=BuildMesh (sphere),Element=Sphere12,Bonds=No DelObj (sphere) Style Ball HideAtom Obj (bud) GlobalZ<50 # Scale the bud down and place it at the cell wall startscale=0.3e0 scale=startscale ScaleObj (bud),X=(scale),Y=(scale),Z=(scale) MoveObj (bud),X=(-scale*radius),Y=(budecc),Z=(-budecc) # Calculate position of bud center in cell coordinates budpos=(cubelen-scale*radius),(cubelen*0.5+budecc),(cubelen*0.5-budecc) # Calculate start volume, the factor (4/3)*Pi is not included since it cancels out startvol=cube (startscale*radius) # Initially everything is fixed FixAll # Initialize simulation ForceField NOVA Longrange None Cutoff 3 TimeStep 2,0.5 SimSteps 10,10 Temp 298K TempCtrl Rescale ChargeAtom all,0 Sim On # Shift the genome in shift=0.05 Console off Pos X=-22,Z=195 Brake 4000 ColorMol Obj 1,blue,cyan ColorAtom PRNA NRNA RNA,Grey HUD Sim # We start by moving all molecules except the growing bud movelist()=ListMol Obj !(bud) #SaveMPG /data/coronamovie2.mp4,1024,720,fps=30,Skip=159,Menu=Yes,JustMacro=Yes lastshot=-1 for i=0. to genlenx step shift Wait 4 MoveMol (join movelist),X=(-shift) if 0+i*0.25!=0+(i+shift)*0.25 # We moved by 4 A, time to adjust the bud # Calculate the new bud scaling factor vol=startvol+(cube radius-startvol)*(i/genlenx) #print 'Vol=(0.00+vol), StartVol=(0.00+startvol), curt=(0.00+curt vol), (0.000+i/genlenx)' newscale=curt vol/radius change=newscale/scale #print 'Scale=(0.000+scale), NewScale=(0.000+newscale), Change=(0.000+change)' # Move the bud center to 0/0/0, scale it, and move it back where it belongs MoveAtom Obj (bud),(-budpos) ScaleObj (bud),X=(change),Y=(change),Z=(change) budpos1=budpos1-(newscale-scale)*radius print (budpos) MoveAtom Obj (bud),(budpos) scale=newscale # Get a new list of automatically moving molecules fixedlist()=ListMol fixed localX>=(cubelen-5) and not Obj (bud) # Free atoms that are new inside the bud FreeMol not (join fixedlist) and not Obj (bud) # Short steepest descent to avoid simulation failures due to overlaps ShowMessage 'Steepest descent minimization at pos (0.00+i)...' TempCtrl SteepDes # Wait until convergence (the maximum speed increases again) lastspeedmax=1e99 while 1 Wait 1 print 'SpeedMax=(speedmax), Last=(lastspeedmax)' if speedmax>lastspeedmax or !speedmax break lastspeedmax=speedmax TempCtrl Rescale # Save single snapshots shot=0+(i*50/genlenx) if shot!=lastshot lastshot=shot SavePNG '/data/coronamovie_(shot).png',Menu=yes freeatoms=CountAtom !fixed density=(0.+freeatoms)/(4.*Pi*cube (scale*radius)/3.) ShowMessage '(100.00*i/genlenx)% completed, (0.000+density)/(SimWarnings)...' Sim Pause #SaveMPG off # EMBED PROTEIN IN MEMBRANE # ========================= # 'solutename' must be present in YOb-format. # 'membranename' with appended 'postfix' is returned by the BuildMem function in md_library.mcr. # The parameters 'cleaned' to 'memextension', 'square', 'confirmneeded' and 'delay' # are explained at the OriInMem function in md_library.mcr. # If 'truncated', the protein will be fixed and temporarily truncated above and below # the membrane to speed up embedding. # 'rhombuslen' is the membrane rhombus side length, or zero if not a rhombus. # The parameters 'ions' and 'cutoff' are explained at the SolvateMem function in md_library.mcr. # If 'shrink', the solute will be shrunken in the end. If 'shrink'>100 the smallest used # pet element will be 'shrink'. # If 'randomized', orientation around the Y-axis and placement inside the membrane # will be randomized. # Membrane will be bent with 'bendradius'. # NOTE: In the end all objects will be removed from the soup and the SimCELL deleted # in preperation of subsequent macro operations, e.g. FillCellList. def EmbedInMem solutename,solutenum,membranename,postfix,cleaned,ph,waterextension,memextension,truncated, square,rhombuslen,confirmneeded,delay,ions,cutoff,shrink,randomized,bendradius global memheight,memcoreheight,PI if !count memheight SetDefaults Remove DelObj SimCELL objname='(solutename)(solutenum)(postfix)' # Check if the solute is already in the soup objs=CountObj '(objname)' if not objs filename='(solutename)(solutenum)_inmem.yob' t1=FileTime (solutename).yob t2=FileTime (filename) if t2>=t1 # Solute has been embedded in membrane before obj=LoadYOb (filename) else # Load and orient solute in membrane obj=LoadYOb (solutename) OriInMem obj,cleaned,ph,waterextension,memextension,square,1,confirmneeded,delay supobj='' if truncated # Truncate protein above and below membrane supobj=DuplicateObj (obj) RemoveObj (supobj) _,memposy=PosObj MembPreview # First select residues in the membrane center and lipid anchors cenlist()=ListRes GlobalY<(memposy+1) and GlobalY>(memposy-1) liplist()=ListRes SmilesMEMBER CH2?CH2CH2CH2CH2? or SmilesMEMBER CH2?CH=C?CH2? and GlobalY<(memposy+memcoreheight*0.5) and GlobalY>(memposy-memcoreheight*0.5) if count cenlist or count liplist reslist()=ListRes (join cenlist) (join liplist) residues=-1 while (residues!=count reslist) # Select iterativly bonded residues inside the membrane core residues=count reslist reslist()=ListRes GlobalY<(memposy+memcoreheight*0.5) and GlobalY>(memposy-memcoreheight*0.5) and with bond to (join reslist) or (join reslist) for i=1 to 5 # Select iterativly bonded residues up to five residues away from residues inside the membrane core # Alpha helix: 5 x 1.5 A per residues=7.5 A=(memheight-memcoreheight)/2 reslist()=ListRes GlobalY<(memposy+memheight*0.5) and GlobalY>(memposy-memheight*0.5) and with bond to (join reslist) or (join reslist) DelRes not (join reslist),Center=Yes Clean PresentObj obj,delay _,memposy=PosObj MembPreview size1,size2,size3,cpos1,cpos2,cpos3=GetMemCell memposy,waterextension,memextension Cell (size) PosObj SimCell,(cpos) # Insert the solute into the membrane randomnum=randomized if solutenum!='' randomnum=randomnum*(solutenum-1) InsertIntoMem obj,membranename,rhombuslen,randomnum,supobj if truncated FixObj (obj) # Solvate system SolvateMem ph,ions,cutoff # Equilibrate membrane TimeStep 2,1.25 Temp 298K Cutoff (cutoff) SimSteps Screen=10,Pairlist=10 EquilibrateMem 5,2.5 Sim Off DelObj Water FreeAll Style Ribbon,Off,Ball LightRes PEA PCH PSE CLR,Y if bendradius # Bend membrane lipids and shift protein down bendsel='Res PEA PCH PSE CLR' Cell Auto,0,Cuboid,not (bendsel) clen1,_,clen2=Cell radius=(max clen)*.5 # Heuristically take only half of calculated shift shift=(cos (radius/bendradius*360/(PI*2))-1)*bendradius*.5 MoveAtom not (bendsel),0,(shift) BendAtom (bendsel),(bendradius) if truncated supobj=LoadYOb (solutename) CleanObj (supobj) SupObj (supobj),(obj),Match=yes SwapObj (supobj),(obj) DelObj (supobj) # Tag membrane as secondary for farcolor ColorFarAtom Obj Membrane,Secondary JoinObj (obj),Membrane,No SwapObj (obj),Membrane if shrink # Additionally shrink the membrane protein petelement=102 if shrink>100 petelement=shrink ShrinkObj (obj),Element=(petelement) ExpandObj (obj) SaveYOb (obj),(filename) NameObj (obj),'(objname)' # Clean up Remove DelObj SimCELL # EMBED A LIST OF PROTEINS INTO A MEMBRANE # ======================================== # Objects listed by name in 'objlist' will be embedded into a lipid bilayer to create # a set of membrane fragments. An object can appear multiple times in the list to create # embeddings with varying orientation and placement in the membrane. These multiple # entries will be returned numbered in 'objlist' to reflect actual object names in the soup. # 'liststep' denotes the number of steps in 'objlist' between two object names, # as 'objlist' can carry additional information about each object, e.g. color. # See function EmbedInMem for a description of remaining parameters. def EmbedListInMem objlist,liststep,membranename,postfix,cleaned,ph,waterextension,memextension,truncated, square,rhombuslen,confirmneeded,delay,ions,cutoff,shrink,bendradius caller (objlist) listlen=count (objlist) for i=1 to listlen step liststep obj=(objlist)(i) num='' randomized=0 oldobjlist(i)=obj cnt=0 for j=1 to i-liststep step liststep if obj==oldobjlist(j) cnt=cnt+1 if cnt # The same fragment is multiple times in the list, enumerate and randomize num=cnt+1 (objlist)(i)='(obj)(num)' randomized=1 EmbedInMem obj,num,membranename,postfix,cleaned,ph,waterextension,memextension,truncated,square,rhombuslen,confirmneeded,delay,ions,cutoff,shrink,randomized,bendradius for i=1 to listlen step liststep (objlist)(i)=(objlist)(i)+postfix # FILL CELL COMPARTMENT WITH LIST OF GROUPS OF PROTEINS # ===================================================== # A cell compartment will be filled with shrunken objects listed in 'objlist'. # 'groupsize' consecutive objects will be grouped together to fill the compartment. # 'liststep' denotes the number of steps in 'objlist' between two object names, # as 'objlist' can carry additional information, e.g. number of copies or color. # Each object color and primary far color will be set with the color in 'objlist'. # The smallest used pet atom element during shrinking is either 'petelement' # or the default pet element if 'petelement'==''. # All parameters for FillCellObj are gathered in a string in 'parameter'. def FillCellGroupList objlist(),groupsize,liststep,petelement,parameter listlen=count objlist groupstep=groupsize*liststep for i=1 to listlen step groupstep copies=0 for j=i to i+groupstep-1 step liststep objname=objlist(j) if strlen '(objname)'>12 # Object names are truncated to 12 characters after loading chrlist()=crack objname objname='' for k=1 to 12 objname=objname+chrlist(k) # Check if object is already in the soup, e.g. after embedding obj=ListObj (objname) if obj atms=CountAtom Obj (obj) if not atms # Add previously removed object to soup AddObj (obj) else # Object not present, load from file obj=LoadYOb (objlist(j)) atms=CountAtom Element 101-112 and Obj (obj) if not atms ShrinkObj (obj),1,(petelement) # By default, fill the compartment with as many copies as possible if liststep>1 if !copies # Each object in 'objlist' is followed by the number of copies copies=objlist(j+1) elif copies!=objlist(j+1) RaiseError 'Object (obj) has not the same number of copies (objlist(j+1)) as object (grouplist1) with (copies) copies, the first member of its group of (groupsize)' if liststep>2 ColorObj (obj),(objlist(j+2)) ColorFarObj (obj),Primary,(objlist(j+2)) grouplist(count grouplist+1)=obj FillCellObj (join grouplist),(copies),(parameter) DelVar grouplist def FillCellList objlist(),liststep,petelement,parameter FillCellGroupList objlist,1,liststep,petelement,parameter # PAIR TWO PROTEINS THAT LINK TWO MEMBRANES # ========================================= # When building a gigastructure that includes two parallel membranes (like the synaptic gap), # there are often transmembrane proteins in both membranes responsible for holding the membranes # together (e.g. Cadherin). Having called EmbedInMem with the protein 'topname' in the top membrane # facing downwards and the protein 'botname' in the bottom membrane facing upwards, you can # dimerize them with this function. 'memdis' is the distance between the membrane meshes that are # filled with proteins, and 'pairlist' is a list of atom pairs and their distance that are known to # interact when the top and bottom proteins dimerize. These will be brought close together by # adding restraints. If 'symmetric', the two proteins are the same, and the restraints in 'pairlist' # are added in both directions. # If 'shrink', the solute will be shrunken in the end. If 'shrink'>100 the smallest used pet element will be 'shrink'. # EXAMPLE for pairlist: 'N Res Asp 1','CD Res Glu 89',3.28, 'N Res Asp 1','OD1 Res Asn 27',2.78,... # NOTE: In the end all objects will be removed from the soup and the SimCELL deleted # in preperation of subsequent macro operations, e.g. FillCellList. def PairInMem topname,botname,memdis,pairlist(),symmetric,shrink global memheight Console off # Remove everything except the two proteins to be paired DelObj SimCell Remove AddObj (topname) (botname) # Make sure that they are present for side in 'top','bot' objlist()=ListObj ((side)name) if count objlist!=1 RaiseError 'Exactly one object named ((side)name) was expected to be present from a previous embedding step, but (count objlist) were found' (side)=objlist(1) # Check the modification times to see if the pairs need to be updated t1=FileTime (topname)_inmem.yob t2=FileTime (topname)_inmempaired.yob t3=FileTime (botname)_inmem.yob t4=FileTime (botname)_inmempaired.yob if t2>t1 and t4>t3 # Pairing already done before, replace with paired objects for side in 'top','bot' DelObj ((side)name) (side)=LoadYOb ((side)name)_inmempaired NameObj ((side)),((side)name) else # Pair the proteins # Convert memdis to float memdis=1.*memdis # Double pairlist if requested if symmetric for i=1 to count pairlist step 3 pairlist(count pairlist+1)=pairlist(i+1) pairlist(count pairlist+1)=pairlist(i) pairlist(count pairlist+1)=pairlist(i+2) # Clear the transformation history PosOriObj (top) (bot),0,0,50,0,0,0 TransformObj (top) (bot),KeepPos=Yes # Move them apart MoveObj (top),Y=(memdis/2) MoveObj (bot),Y=(-memdis/2) # Do they clash? This could cause unresolvable knots during energy minimization bumps=CountAtom Obj (top) with distance<2 from Obj (bot) if bumps>3 RaiseError 'Unfortunately (topname) and (botname) bump into each other initially, please change one conformation manually' # Prepare simulation Clean Cell Auto,Extension=10 ZoomObj SimCell,Steps=1 ForceField NOVA,SetPar=Yes SimSteps 10,10 # Fix the membrane parts FixAtom GlobalY<((-memdis+memheight)/2) or GlobalY>((memdis-memheight)/2) # Constrain all hydrogen bonds to stabilize the model if Structure OptHydAll ShowMessage 'Stabilizing secondary structure with distance restraints...' Wait 1 for side in 'top','bot' hbolist()=ListHBoAtom Obj ((side)),Obj ((side)),Results=6 for i=1 to count hbolist step 6 # Add a spring between the acceptor (either hbolist(i) or hbolist(i+1)) and the hydrogen (hbolist(i+5)) AddSpring (hbolist(i+(hbolist(i+3)=='don'))),(hbolist(i+4)),Len=(hbolist(i+5)),SFC=70 # Constrain the distances with increasing strength (stretching force constants in N/m) sfc=0. while sfc<50 # Add user-provided distance restraints for j=1 to count pairlist step 3 #print (pairlist(j)) (pairlist(j+1)) (pairlist(j+2)) topsel=pairlist(j) botsel=pairlist(j+1) # Get the atom numbers and verify that each selection matches exactly 1 for side in 'top','bot' (side)list()=ListAtom ((side)sel) Obj ((side)) if count (side)list!=1 RaiseError 'Found (count (side)list) instead of one expected atom ((side)sel) in object ((side)name)' # Add the distance restraint AddSpring (toplist1),(botlist1),Len=(pairlist(j+2)),SFC=(sfc) # Do a steepest descent minimization ShowMessage 'Steepest descent minimization with restraint force constant (sfc)' TempCtrl SteepDes Sim On Wait 10 for j=1 to 100 Wait 1 if SpeedMax<4000 break # And some MD ShowMessage 'Simulated annealing minimization with restraint force constant (sfc)' TempCtrl Anneal Temp 298K Wait 10 for j=1 to 100 Wait 1 if SpeedMax<500 break # Increase SFC if !sfc sfc=1e-2 else sfc=sfc*1.5 # Save paired structures and undo the transformations applied during the simulation Free for side in 'top','bot' SaveYOb ((side)),((side)name)_inmempaired,Transform=Yes if shrink # Additionally shrink the membrane protein DelObj ((side)) (side)=LoadYOb ((side)name)_inmempaired NameObj ((side)),((side)name) petelement=102 if shrink>100 petelement=shrink ShrinkObj ((side)),Element=(petelement) ExpandObj ((side)) SaveYOb ((side)),((side)name)_inmempaired # Clean up Remove DelObj SimCELL # LOAD OBJECT POSITIONS AND ORIENTATIONS FROM FILE # ================================================ def LoadPacking filename datatab=LoadTab (filename) # Count the total number of proteins proteinlist()=Tab (datatab),Column=1 proteins=count proteinlist for i=1 to proteins # For every line in the file get the corresponding transform data and apply them to the instanced protein protein,posori()=Tab (datatab),Row=(i) obj=ListObj (protein) if obj # Object is already in the soup atms=CountAtom Obj (obj) if not atms # Add previously removed object to soup AddObj (obj) # Add instance InstanceObj (obj),copies=1 else obj=LoadYOb (protein) ins=InstanceObj (obj) # Place instances and convert from rightHanded xzy to leftHanded xzy PosOriObj (obj) Ins (ins),(posori1),(posori2),(-posori3),(posori4),(180.-posori5),(-posori6) if not (i%1000) ShowMessage 'Packing (protein), scene (i*100./proteins)% completed...' Wait 1 HideMessage