read_swc (v0.9.0, graph backend) gives 2 compartments for soma despite ncomp=1

[kirkw-1]

The attached spreadsheet contains first few lines of an SWC morphology, the corresponding Jaxley constructed morphology, and the read_swc() code used to import. My question is, does the soma construction make sense according to your understanding of how the import should work?
My ncomp is 1, but jaxley assigned my soma 2 compartments. I thought this was related to the -1 root node issue, but read in the github that version 0.9.0 removed that step. And if it were related to root node, why is the length 10micron?
Thanks,
jaxley output morphology.ods

[michaeldeistler]

Hi Kirk!

Thanks for creating this!

Two things:

1. Yes, it makes sense to me that there are two somatic compartments. The reason is that the axon is connected to the middle of the somatic traced points. As such, Jaxley splits the soma into two parts (separated by the branchpoint to the axon). This is not related to the root node. Does this make sense to you?

2. In v0.9.0, I get different a length than the length that you report (but still two somatic compartments, see below). In particular, I get:

   global_comp_index     length    radius   axon  basal   soma
0                  0  10.426234  6.502691  False  False   True
1                  1  12.743175  6.728561  False  False   True
2                  2  60.000000  1.000000   True  False  False
3                  3   8.254651  0.558292  False   True  False

I could imagine that this is related to this issue, which reported issues when upgrading to v0.9.0. In that issue, the solution was to start a new conda environment from scratch and reinstall Jaxley. I have to admit that I am not sure what is causing this issue. You suggest you first update the following packages:

pip install pandas --upgrade
pip install numpy --upgrade
pip install jax --upgrade
pip install jaxlib --upgrade 
pip install networkx --upgrade

If this does not help, then I recommend starting a new conda environment.

Hope this helps!
Michael

[kirkw-1]

did you mean import CaNernstReversal from l5pc.py? I have latest jaxley_mech 0.3.0 and can't find any CaNernstReversal file. If so, then all my simulations have included that import l5pc.py so should have CaNernstReversal already, and now I just tested a toy model based on your 01_synthetic from the paper, and it also has cell.insert(CaNernstReversal()) but also returns eca 0.

[michaeldeistler]

I unfortunately cannot reproduce the issue (but yes, I meant importing from l5pc.py; I fixed it in my message above). Can you try:

from jax import config

config.update("jax_enable_x64", True)
config.update("jax_platform_name", "cpu")

import os
os.environ["XLA_PYTHON_CLIENT_MEM_FRACTION"] = ".8"

import matplotlib.pyplot as plt
import numpy as np
import jax.numpy as jnp

import jaxley as jx
from jaxley.channels import Leak
from jaxley_mech.channels.l5pc import *

gt_apical = {}
gt_soma = {}
gt_axon = {}
gt_apical["apical_NaTs2T_gNaTs2T"] = 0.026145
gt_apical["apical_SKv3_1_gSKv3_1"] = 0.004226
gt_apical["apical_M_gM"] = 0.000143
gt_soma["somatic_NaTs2T_gNaTs2T"] = 0.983955
gt_soma["somatic_SKv3_1_gSKv3_1"] = 0.30347
gt_soma["somatic_SKE2_gSKE2"] = 0.008407
gt_soma["somatic_CaPump_gamma"] = 0.000609
gt_soma["somatic_CaPump_decay"] = 210.485291
gt_soma["somatic_CaHVA_gCaHVA"] = 0.000994
gt_soma["somatic_CaLVA_gCaLVA"] = 0.000333
gt_axon["axonal_NaTaT_gNaTaT"] = 3.137968
gt_axon["axonal_KPst_gKPst"] = 0.973538
gt_axon["axonal_KTst_gKTst"] = 0.089259
gt_axon["axonal_SKE2_gSKE2"] = 0.007104
gt_axon["axonal_SKv3_1_gSKv3_1"] = 1.021945
gt_axon["axonal_CaHVA_gCaHVA"] = 0.00099
gt_axon["axonal_CaLVA_gCaLVA"] = 0.008752
gt_axon["axonal_CaPump_gamma"] = 0.00291
gt_axon["axonal_CaPump_decay"] = 287.19873

comp = jx.Compartment()
branch = jx.Branch(comp, 4)
cell = jx.Cell(branch, [-1, 0, 0, 0])
cell.branch(0).add_to_group("soma")
cell.branch(1).add_to_group("basal")
cell.branch(2).add_to_group("apical")
cell.branch(3).add_to_group("axon")

########## APICAL ##########
cell.apical.insert(NaTs2T().change_name("apical_NaTs2T"))
cell.apical.insert(SKv3_1().change_name("apical_SKv3_1"))
cell.apical.insert(M().change_name("apical_M"))
cell.apical.insert(H().change_name("apical_H"))

########## SOMA ##########
cell.soma.insert(NaTs2T().change_name("somatic_NaTs2T"))
cell.soma.insert(SKv3_1().change_name("somatic_SKv3_1"))
cell.soma.insert(SKE2().change_name("somatic_SKE2"))
ca_dynamics = CaNernstReversal()
ca_dynamics.channel_constants["T"] = 307.15
cell.soma.insert(ca_dynamics)
cell.soma.insert(CaPump().change_name("somatic_CaPump"))
cell.soma.insert(CaHVA().change_name("somatic_CaHVA"))
cell.soma.insert(CaLVA().change_name("somatic_CaLVA"))
cell.soma.set("CaCon_i", 5e-05)
cell.soma.set("CaCon_e", 2.0)

########## BASAL ##########
cell.basal.insert(H().change_name("basal_H"))
cell.basal.set("basal_H_gH", 8e-5)

# ########## AXON ##########
cell.insert(CaNernstReversal())
cell.set("CaCon_i", 5e-05)
cell.set("CaCon_e", 2.0)

cell.axon.insert(NaTaT().change_name("axonal_NaTaT"))
cell.axon.set("axonal_NaTaT_gNaTaT", 3.137968)

cell.axon.insert(KTst().change_name("axonal_KTst"))
cell.axon.set("axonal_KTst_gKTst", 0.089259)

cell.axon.insert(CaPump().change_name("axonal_CaPump"))
cell.axon.set("axonal_CaPump_gamma", 0.00291)
cell.axon.set("axonal_CaPump_decay", 287.19873)

cell.axon.insert(SKE2().change_name("axonal_SKE2"))
cell.axon.set("axonal_SKE2_gSKE2", 0.007104)

cell.axon.insert(CaHVA().change_name("axonal_CaHVA"))
cell.axon.set("axonal_CaHVA_gCaHVA", 0.00099)

cell.axon.insert(KPst().change_name("axonal_KPst"))
cell.axon.set("axonal_KPst_gKPst", 0.973538)

cell.axon.insert(SKv3_1().change_name("axonal_SKv3_1"))
cell.axon.set("axonal_SKv3_1_gSKv3_1", 1.021945)

cell.axon.insert(CaLVA().change_name("axonal_CaLVA"))
cell.axon.set("axonal_CaLVA_gCaLVA", 0.008752)


########## WHOLE CELL  ##########
cell.insert(Leak())
cell.set("Leak_gLeak", 3e-05)
cell.set("Leak_eLeak", -75.0)

cell.set("axial_resistivity", 100.0)
cell.set("eNa", 50.0)
cell.set("eK", -85.0)
cell.set("v", -65.0)


for key in gt_apical.keys():
    cell.apical.set(key, gt_apical[key])

for key in gt_soma.keys():
    cell.soma.set(key, gt_soma[key])

for key in gt_axon.keys():
    cell.axon.set(key, gt_axon[key])

dt = 0.025
t_max = 100.0
time_vec = np.arange(0, t_max+2*dt, dt)

cell.delete_stimuli()
cell.delete_recordings()

i_delay = 5.0
i_dur = 90.0
i_amp = 3.0
current = jx.step_current(i_delay, i_dur, i_amp, dt, t_max)
cell.soma.branch(0).loc(0.0).stimulate(current)
cell.soma.branch(0).loc(0.0).record("eCa")

cell.set("v", -65.0)
cell.init_states()

eca = jx.integrate(cell)
plt.plot(time_vec, eca[0])

For me, this plots the following:

Please let me know if you get the same results with this code!

Michael

[kirkw-1]

Thanks, yes I just ran your code and it worked fine, and I've confirmed that eca is definitely working ok in my model. But I still get eca 0.0 when I run cell.show() after the simulation, is that expected?

[michaeldeistler]

Yes, that is expected. .show() does not update the states after a simulation---it shows the values that are used to start a simulation. If you want to get all final states (after the simulation), do:

recordings, final_states = jx.integrate(cell, return_states=True)

This returns the final_states as a Dictionary. This will still not be shown in .show() though. For them to be shown in .show(), you have to do the following:

cell.set("eCa", final_states["eCa"][cell._internal_node_inds])

This is a bit clunky, and I opened an issue here to improve this.

Hope this helps!
Michael

[kirkw-1]

ah I see, OK thanks very much for your help, that explains it.