Counting T1 transitions

[glyg]

Here is an email from Chi Xu, who works with @gcourcou

This is Chi, I am working with George Courcoubetis on tissue research. Currently we want to visualize T1 transitions in a growing tissue. Specifically, we want to get the number of T1 transitions and the place(which cells) where they happen in the tissue when we use find energy min function with T1 transitions is on. However I found it's not easy to do so. I keep track of the number of sides for each cell, if it's changing I assume T1 transitions happen. But it could be different(like cell divisions, or T1 transitions happen twice on a cell which keeps the number of sides the same for that cell). So I am sending this email to ask whether there is an explicit way to show the number of T1 transitions and the place where they happen during finding energy min process.

For now there is no explicit way to do that.

A solution would be to be able to pass a on_topo_change callback to the QSSolver class to be called here. I should have time to open a PR for that next week. It can be done by subclasseing QSSolver:

class QSSolverT(QSSolver):

    def __init__(self, on_topo_change, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.on_topo_change = on_topo_change

    def _minimize(self, eptm, geom, model, **kwargs):
        for i in count():
            if i == MAX_ITER:
                return self.res
            pos0 = eptm.vert_df.loc[
                eptm.vert_df.is_active.astype(bool), eptm.coords
            ].values.flatten()
            try:
                self.res = optimize.minimize(
                    self._opt_energy,
                    pos0,
                    args=(eptm, geom, model),
                    jac=self._opt_grad,
                    **kwargs
                )
                return self.res
            except TopologyChangeError:
                log.info("TopologyChange")
                self.on_topo_change(eptm)
                self.num_restarts = i + 1

A longer term solution, also more generic & interesting is to have a better History object that could keep track of the objects throughout rearrangements (e.g as a 4D graph). This could go with a redefinition of the I/O.

[xcricardo]

Thank you Guillaume! What is on_topo_change? How can we get the number of T1 transisions by using it? Is the code above enough to obtain the number of T1 transitions?

[glyg]

Hi, on_topo_change should be a function taking eptm as it's only parameter. For example:

def count_t1(eptm):
    if "num_type1" in eptm.settings:
        eptm.settings["num_type1"] += 1
    else:
        eptm.settings["num_type1"] = 1

solver = QSSolverT(on_topo_change=count_t1)
. . .

I realize now that this would count more changes than just T1s...

Maybe a better way to do it would be through the "opposite_face" column in sheet.face_df. The get_opposite_face method would allow you to see if there are neighbor changes.

[gcourcou]

Regarding the T1 transitions: We wrote a code that outputs T1 transitions. We did not identify any T1 transitions during our growing tissues. We use the auto T1 transition option in the energy minimization step. Have you tested that auto_t1 actually works? Also we noticed a magic number to identify short edges eligible for T1 transitions 1e-6. Why was that choice made? Can we vary it and see its effect? Thank you so much, we currently submitted to PLOS Comp bio and we have positive reviews but need to address many questions. It is a bit worrisome that we do not see any T1s, it would seem like a great way for the tissue to relax, which could lead to the prevention of avalanches, which is the central subject of our manuscript.

[xcricardo]

This is our code for testing the t1 transitions between frames. It failed when the indexes got shuffled.

for jump_index in jump_indexes:
    print(jump_index)
    dsets1 = hdf5.load_datasets(name + str(jump_index - 1) + ".hf5")
    specs1 = config.geometry.planar_sheet()
    before = Sheet("periodic", dsets1, specs1)

    dsets2 = hdf5.load_datasets(name + str(jump_index) + ".hf5")
    specs2 = config.geometry.planar_sheet()
    after = Sheet("periodic", dsets2, specs2)

    # t1 transitions
    # find neighborlist first
    beforeneighbor = dict()
    afterneighbor = dict()
    for index, row in before.edge_df.iterrows():
        key = row["face"]
        val = row["opposite"]
        if val != -1:
            if key in beforeneighbor:
                beforeneighbor[key].add(before.edge_df.at[val, "face"])
            else:
                beforeneighbor[key] = set()
                beforeneighbor[key].add(before.edge_df.at[val, "face"])

    for index, row in after.edge_df.iterrows():
        key = row["face"]
        val = row["opposite"]
        if val != -1:
            if key in afterneighbor:
                afterneighbor[key].add(after.edge_df.at[val, "face"])
            else:
                afterneighbor[key] = set()
                afterneighbor[key].add(after.edge_df.at[val, "face"])

    # compare
    #    difface = []
    #    for index, row in before.face_df.iterrows():
    #        if row['num_sides'] != after.face_df.at[index, 'num_sides']:
    #            difface.append(index)
    #    print("difface is:")
    #    print(difface)
    t1_set = set()
    for key in afterneighbor:
        if key in beforeneighbor:
            for face in afterneighbor[key]:
                if face not in beforeneighbor[key]:
                    #                    print("new neighbor:")
                    #                    print(face)
                    if face in beforeneighbor:
                        print("t1:")
                        print(face)
                        t1_set.add(key)`

[glyg]

@gcourcou :

Regarding the T1 transitions: We wrote a code that outputs T1 transitions. We did not identify any T1 transitions during our growing tissues. We use the auto T1 transition option in the energy minimization step. Have you tested that auto_t1 actually works?

Reading the code, it should work but I haven't tested it in some time. I'll draft a test today to ensure T1 do happen. I checked, and they do, see link to notebook at the end of this answer

Also we noticed a magic number to identify short edges eligible for T1 transitions 1e-6.

Yes you can see it here

Why was that choice made?

This is a default value, sufficiently small to avoid problems in case of small tissue. Other than that it is arbitrary but

Can we vary it and see its effect?

Absolutely, it is the "threshold_length" parameter, its actual value is read from the sheet.settings dictionnary (as you can see in the code line linked above). So setting:

sheet.settings['threshold_length'] = 1e-2

will change the threshold fot T1 transitions.

, we currently submitted to PLOS Comp bio and we have positive reviews

Awesome, congrats!

It is a bit worrisome that we do not see any T1s, it would seem like a great way for the tissue to relax, which could lead to the prevention of avalanches, which is the central subject of our manuscript.

Hopefully changing the threshold will do the tick, I'll try to have a test ready quickly

---

@xcricardo

The spirit is good, should not take much to fix. As you noticed, re-indexing scrambles the indices, so in order to have consistent data, you need to add an "id" column to the objects you want to track:

sheet.face_df['id'] = sheet.face_df.index

(This column is later updated if present when a cell division occurs)

If you don't have division or apoptosis though, face indices should not change (I don't know if it's your case).

I cooked up a function that returns the set of face pairs that changed between to transitions:
Look at this notebook

Hope all of this helps!

[gcourcou]

Thank you Guillaume yes, this helps!! We do have mitosis 100% We will look into your notebook and see if we can use parts of it.

[glyg]

So yes if you have mitoses, you need to setup the "id" at the bigining of the simulation. This actually should be the default for cells (faces in 2D) I think, it never makes sens to not want that info!

I think with the linked code, you can later on make sense of wether a T1 happened or a division

[gcourcou]

Thanks!! I do want to avoid re-running the simulations but it should be easy to do!

[gcourcou]

@glyg I used

 sheet.face_df['id'] = sheet.face_df.index

and I am getting duplicate IDs in my sheet.face_df['id'] column. E.g. there are multiple cells with an id of 2 e.t.c. I use daughter = cell_division(sheet, index, geom, angle=angle_div) for the division. (from tyssue.topology.sheet_topology import cell_division)

Let me know how I can fix that. Thanks again for all your help!!

[glyg]

Ha I thought the id column was updated in cell_division, but it must not be the case, sorry. So you have to manually do:

daughter = cell_division(sheet, index, geom, angle=angle_div)
sheet.face_df.loc[daughter, "id"] = sheet.face_df.index.max() + 1

[glyg]

Also for now as a bug fix (in case you don't want to re-run your simulations!) , daughter cells will have the same id as their mother but I think always a higher index (re-indexing only shifts indices to the right when a cell is missing for example), so you should be able to de-uniquify them.

[gcourcou]

Thanks @glyg , really appreciate your prompt responses! Lastly, could you elaborate on how auto_t1 works in relation to energy and IH transition (which i do not understand). Is it basically saying if an edge is less than the threshold, then try to do a T1, if the T1 leads to lower energy, then keep it?

[glyg]

Hi, no problem,

No, the T1 just happens if the length is lower than the threshold, there is no energy evaluation. The gradient descent is simply restarted. The transition would be reverted only if the length of the new junction stays lower at the next gradient descent step.

This is actually a problem with T1 transitions described like that, the possibility of oscillation. The Finengan et al. paper has a better algorithm (which is implemented in tyssue, but only integrated in the EulerSolver, it is less straight-forward to do so without a time scale in a quasistatic model): T1 is not a single step event but rather first you form a 4-way vertex, that gets resolved stochastically by one of the 4 junctions detaching.

You can look at this notebook for more details

[gcourcou]

Ha I thought the id column was updated in cell_division, but it must not be the case, sorry. So you have to manually do:

daughter = cell_division(sheet, index, geom, angle=angle_div)
sheet.face_df.loc[daughter, "id"] = sheet.face_df.index.max() + 1

We are still getting repeating indexes in the output. Any guess why ?
Maybe use sheet.face_df.id.max() instead of .index.max()?

[glyg]

Maybe use sheet.face_df.id.max() instead of .index.max()?

That might be a good idea yes :/

[gcourcou]

Hi @glyg could you explain how auto_t3 works in the code? (should be t2 not t3! according to literature) What is the parameter that sets if a cell will be removed?

[glyg]

Hi @gcourcou sure I can explain. relevant code is
here

We follow Okuda et al (the Recombinant Network Reconnection paper) rule for a type 3 transition i.e the elimination of a triangular cell from the apical surface -- maybe it's type 2, I haven't revised the literature in a while, and I might have been wrong when I wrote that!
The rule states that a triangular face is eliminated if it's longest edge is shorter than the threshold length, the same setting as for type 1.

[glyg]

and yes according to that paper it is a type 2 transition. We should duplicate the function with a deprecation waring for the old name - implementing the true T3 could be interesting too! Thanks for pointing this out

[gcourcou]

Thanks @glyg !! I appreciate your response. It is slightly different from the area threshold used in other studies but seems equally valid to me. Glad i could point out something to improve the library. I also would recommend setting the T1 transition length to 0.01 but that is up to you see this link .

[glyg]

Thanks for the ref, I guess we can trust A. Fletcher ^^'

If you don't mind, could you change the default for threshold_length in your opened PR?

[gcourcou]

Sure, done!