Create multi-strand system for mrDNA simulations

[FionaSander]

The final aim is to have a system of several strands which are not connected with each other. Each strand is made up of several doublestranded and singlestranded sections (script)

1. A system of one strand of dsDNA and ssDNA sections is created with drawNA

2. The drawNA system is written into oxDNA files (.conf, .top)

3. The oxDNA files are converted to .pdb file format for use in mrDNA simulation

4. The strand is used for a short coarse-grained simulation in mrdna.

5. The .pdb output file from this initial simulation is converted back into oxDNA files and read back into drawNA

6. The strand is then added to the mother system (this will be the multi-strand system)

7. The simulated strand is copied and shifted in space, converted to oxDNA files and .pdb file format to be read into mrDNA for a short simulation

steps 4 -7 are repeated until x number of strands are described in the mother system.

8. The mother system is the final system of interest and will be simulated with mrDNA (conversion to oxDNA and .pdb file format first).

ISSUE: Addition of strands to mother system. It looks like the strands are not added properly to the system (even just one strand). There might be the need for an additional script for addition of strands into one new mother system.

One strand system without reading into a mother system (works well) - (displayed with ovito and vmd):

One strand system with reading into mother system as a starting point for a multi-strand system - (displayed with ovito and vmd):

[debeshmandal]

@FionaSander I have a question:

Is there a chance the the strand leaves the periodic box during the simulation?

[FionaSander]

@debeshmandal - It definitely looks like there might be some issue with the box size. However, when I run the same system of one strand with the working script I don't get any box size issues. I compared the .conf files of the standard system with the mother system and they are not equal even though they should be when I'm only trying to insert one strand into the mother system.

Do you think this has something to do with how I import the .pdb output from mrdna back into drawNA?

[debeshmandal]

In lines 149-151

#Convert oxDNA output files to .pdb input file format to run mrdna on them
subprocess.call(
     "bash -c \"source ~/.bashrc && module load CUDA && cd {} && python3 /home/fiona/tacoxDNA-python3/src/oxDNA_PDB.py oxdna.turns.top oxdna.turns.conf 35\"".format(file_path), shell=True)

Check that the PDB outputted is in the 3'-5' direction - I think that it is more likely the other way around

[debeshmandal]

If the direction is wrong, it will completed break your system because the nucleotides will want to completely reorient themselves

[debeshmandal]

Let me try and write some pseudo code to track the best algorithm for doing this:

# Fiona Code

DELTA = np.array([0., 5., 0.])
TEMPLATE = 'main.{}'
N = 10
SIMULATIONS = range(1, N)

predefined_functions:
    generate_system
    convert_to_pdb
    run_simulation
    read_pdb
    translate_system

clone_system = generate_system()
clone_name = TEMPLATE.format(0)
clone_pdb = f'{clone_name}.pdb'
system.write_oxDNA(clone_name)
convert_to_pdb(system, fname=clone_pdb)

for i in SIMULATIONS:

    new_name = TEMPLATE.format(i)
    new_pdb = f'{new_name}.pdb'

    run_simulation(clone_pdb, new_pdb)

    # use predefined new_* to define filenames and then import
    # this overwrites the previously defined clone_* variables
    clone_name = new_name
    clone_pdb = new_pdb
    clone_system = read_pdb(clone_pdb)

    # do things with the clone system
    translate_system(clone_system, DELTA)

    # standard operations  to export everything correctly
    # and run simulation
    clone_system.write_oxDNA(clone_name)
    convert_to_pdb(clone_system, fname=clone_pdb)

The result of this is that you should have a load of names that look like this:

oxdna.main.0.conf
oxdna.main.0.top
main.0.pdb
oxdna.main.1.conf
oxdna.main.1.top
main.1.pdb
oxdna.main.2.conf
oxdna.main.2.top
main.2.pdb

And all that is required is that you need to set the TEMPLATE at the beginning

The constants at the beginning should be set in the ArgumentParser object

[debeshmandal]

For a more detailed look at how to set the template:

from pathlib import Path
target_directory = Path('/home/fiona/simulations/target/')
TEMPLATE = str(target_directory / f"l{length}-ds{stapled}-p{percentage}.{{}}")

I suggest keeping this script as one file for one configuration - this means that you can run each configuration as a separate job and their time will be limited.

[FionaSander]

The PDB output is in 3'-5' direction so I think the problem must be somewhere else. I just ran a short simulation giving me .pdb, .conf and .top files for the standard system ("turns") as well as for the mother system ("mother_turns"). The only difference I can find is in the configuration file regarding the box size but when trying to force another box size the strand is still messed up.

Regarding the code above: This clones and translates the strand but does it also add it to a overall system which in the end would contain all modified strands?

mrdna_output.zip

[debeshmandal]

Ah sorry, I did not realise you wanted to copy the system into itself - I will amend the code

[debeshmandal]

@FionaSander I'm curious as to why you are simulating everything in the same simulation? is it only for viewing or are you simulating different stages in the same box

If the latter, then you may run into problems when the strands start to entangle during the simulation due to the thermostat

[debeshmandal]

For having a master system used for viewing only:

# Fiona Code

DELTA = np.array([0., 5., 0.])
MAIN = 'main'
TEMPLATE = MAIN + '.{}'
N = 10
SIMULATIONS = range(1, N)

predefined_functions:
    generate_system
    convert_to_pdb
    run_simulation
    read_pdb
    translate_system


main_system = generate_system()

clone_system = main_system.copy()
clone_name = TEMPLATE.format(0)
clone_pdb = f'{clone_name}.pdb'


system.write_oxDNA(clone_name)
convert_to_pdb(system, fname=clone_pdb)

for i in SIMULATIONS:

    new_name = TEMPLATE.format(i)
    new_pdb = f'{new_name}.pdb'

    run_simulation(clone_pdb, new_pdb)

    # use predefined new_* to define filenames and then import
    # this overwrites the previously defined clone_* variables
    clone_name = new_name
    clone_pdb = new_pdb
    clone_system = read_pdb(clone_pdb)

    # do things with the clone system
    translate_system(clone_system, DELTA)

    # clone the strands and add them to the main system
    clone_strands = clone_system.strands
    main_system.add_strands(clone_strands)
    
    # standard operations  to export everything correctly
    # and run simulation
    clone_system.write_oxDNA(clone_name)
    convert_to_pdb(clone_system, fname=clone_pdb)

main_system.write_oxDNA(MAIN)

[FionaSander]

Ah okay, I thought I missed something! I'll have a look at the code above.

I want to have all the strands in one system because that's what I'm trying to simulate. Similar to the original idea with the polymers cross-linked with the liquid crystals there should be many DNA strands in this system that would naturally entangle etc. over time. The idea is to be able to modify the strand density to achieve alignment of the dsDNA section of my system.

[debeshmandal]

Okay have a look at the above and get back to me on here about any further amendments

[FionaSander]

Will do that, thanks!

[debeshmandal]

In my last message with a full script, files called oxdna.main.conf and oxdna.main.top will be outputted which could be the start of a system that you need to simulate everything together

Also, if there are anythings missing from the code that you think would be helpful, make a note here and I will try to add them

[FionaSander]

@debeshmandal I integrated your suggestions into my code and I can now add different strands to the main system. However, there seems to be an issue with how the double stranded sections are interpreted (see attached screenshot). I have the suspicion it has something to do with the translation using pos_com (translate_system - see below) or how the .pdb file is read into drawNA (read_pdb).
Have you seen anything similar in one of your systems?

    for nt in system.nucleotides:
        nt.pos_com += translation_vector

    return system

I will keep trying to find the source of the detachment of the double stranded section from the ssDNA but let me know if you have any idea on this.

See updated code (lines 150-250).

[debeshmandal]

Try changing this:

# clone the strands and add them to the main system
    clone_strands = clone_system.strands
    main_system.add_strands(clone_strands)

to this:

# clone the strands and add them to the main system
    clone_strands = [i.copy() for i in clone_system.strands]
    main_system.add_strands(clone_strands)

[debeshmandal]

If that doesn't work then I will make a full script and put it in the protocols folder where a short simulation is done using oxDNA only but the translations would be the same, and maybe that will help

[FionaSander]

I tried the above and it doesn't seem to make any difference. I went into the details of the code again and I'm pretty sure the problem must be the translation of the nucleotides as I cannot find any errors in the generated pdb files.

Would you be able to briefly discuss the translation of bases on Teams one of the next days?

[debeshmandal]

I am about to merge some changes to the master branch which should allow transformations to occur. To do translations:

vector = np.array([10., 0., 0.])
nucleotide.translate(vector)
strand.translate(vector)
system.translate(vector)

All of the above will perform translations on the positions of their constituent objects

[FionaSander]

Perfect, I just had a quick look and ran some simulations.
The translate function is still basically the same than what I already implemented so I assume you suggest to use the combination of translate, transform and rotate in test_nucleotide.py/ test_strand.py/ test_system.py ?

[debeshmandal]

I am updating a final protocol which shows how you can simply use system.copy() and system.translate(vector) to create the replicas - trying running the protocol which is found at drawNA/protocols/system-translations to check that you understand how it works (if it's not there, it will be soon, i am just updating'

[debeshmandal]

I think that this should be completed now, so let me know if I can close the issue

[FionaSander]

I'm having a look at this right now.
What is the OXDNA parameter used for when running main.py and calling run_simulation? In both you also give the input file(s) so I'm not sure what input is needed here? Is this just for when one wants to use oxDNA for the simulations? As you say in the notes - I'll adjust the run_simulation function for mrdna simulations and if the oxDNA parameter is not really needed for this I'll just leave it away.

[debeshmandal]

I'm just adding a better description but the oxDNA part refers to the path to the oxDNA binary file and the input is for an oxDNA input file - I know this isn't related to mrDNA but it demonstrates how to do the translation and make replicas

[debeshmandal]

Closing this for now - reopen if any more help is needed

[debeshmandal]

@FionaSander You mentioned in an e-mail that you get the following error:

File "/home/fiona/miniconda3/lib/python3.7/site-packages/mrdna-0.0.0-py3.7.egg/mrdna/segmentmodel.py", line 1226, in _generate_one_bead
    assert(np.linalg.norm(opos-pos) < 10 )
AssertionError

I'm going to try to find the answers to these questions:

1. what is opos?
2. what is pos?
3. why has a value of 10 been used for the condition?

[debeshmandal]

For reference, here is the function

class DoubleStrandedSegment(Segment):
    ...
    def _generate_one_bead(self, contour_position, nts):
        pos = self.contour_to_position(contour_position)
        if self.local_twist:
            orientation = self.contour_to_orientation(contour_position)
            if orientation is None:
                print("WARNING: local_twist is True, but orientation is None; using identity")
                orientation = np.eye(3)
            opos = pos + orientation.dot( np.array((Segment.orientation_bond.r0,0,0)) )
            # if np.linalg.norm(pos) > 1e3:
            #     pdb.set_trace()
            assert(np.linalg.norm(opos-pos) < 10 )
            o = SegmentParticle( Segment.orientation_particle, opos, name="O",
                                 contour_position = contour_position,
                                 num_nt=nts, parent=self )
            bead = SegmentParticle( Segment.dsDNA_particle, pos, name="DNA",
                                    num_nt=nts, parent=self, 
                                    orientation_bead=o,
                                    contour_position=contour_position )

        else:
            bead = SegmentParticle( Segment.dsDNA_particle, pos, name="DNA",
                                    num_nt=nts, parent=self,
                                    contour_position=contour_position )
        self._add_bead(bead)
        return bead

[debeshmandal]

There is an assertion to ensure that opos and pos are within 10 units of each other.

• pos is the position of the contour?
• opos is the position of the orientation particle`

The way that opos is determined is by adding pos to a new vector called orientation.dot(np.array((Segment.orientation_bond.r0, 0, 0) where:

• orientation is a 3x3 matrix
• the argument to npndarray.dot() is array(Segment.orientation_bond.r0, 0, 0)

Segment.orientation_bond.r0 is given by this:

orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )

Therefore one would assume that the equilibrium length of this bond is 1.5

This means that the value of the matrix function looks something like this:

orientation.dot(np.array([1,5, 0, 0]))

Since the norm is greater than 10 there is clearly an issue with this object named orientation:

orientation = self.contour_to_orientation(contour_position)
class Segment:
    ...
    def contour_to_orientation(self,s):
        assert( isinstance(s,float) or isinstance(s,int) or len(s) == 1 )   # TODO make vectorized version

        if self.quaternion_spline_params is None:
            axis = self.contour_to_tangent(s)
            axis = axis / np.linalg.norm(axis)
            rotAxis = np.cross(axis,np.array((0,0,1)))
            rotAxisL = np.linalg.norm(rotAxis)
            zAxis = np.array((0,0,1))

            if rotAxisL > 0.001:
                theta = np.arcsin(rotAxisL) * 180/np.pi
                if axis.dot(zAxis) < 0: theta = 180-theta
                orientation0 = rotationAboutAxis( rotAxis/rotAxisL, theta, normalizeAxis=False ).T
            else:
                orientation0 = np.eye(3) if axis.dot(zAxis) > 0 else \
                               rotationAboutAxis( np.array((1,0,0)), 180, normalizeAxis=False )
            if self.start_orientation is not None:
                orientation0 = orientation0.dot(self.start_orientation)

            orientation = rotationAboutAxis( axis, self.twist_per_nt*self.contour_to_nt_pos(s), normalizeAxis=False )
            orientation = orientation.dot(orientation0)
        else:
            q = interpolate.splev( s, self.quaternion_spline_params[0] )
            if len(q) > 1: q = np.array(q).T # TODO: is this needed?
            orientation = quaternion_to_matrix(q)

        return orientation

[debeshmandal]

Clearly the orientation should be a matrix that has an np.linalg.norm of 1.0 but this isn't happening, why this the case I have no idea because in our code, we check to make sure that this is the case - and should not be caused by an issue with drawNA

However, I did find an issue when creating replicas that sometimes they were outside of the box, I think that it would be useful to check that everything is in the box, because if it is not, when a new particle is created, it could be moved to the opposite side of the box because of periodic boundary conditions, hence the distance between new (created) particle and the old (imported) particle would be the (length of the box) - (ideal bond length)

Overall, I think that this is a problem with the box more than anything else - and you didn't see this problem before because the double-stranded parts of the system were not being copied properly

Let me know how editing the box goes on here please