Modeling the multiplication of two binary decision variables in undirected graph in python

In an undirected graph, I’m trying to model a constraint that forcing the optimizer to set an edge $(u,v)$ between two nodes to only exist (= $1$) if the two nodes have been selected to be $1$. The three decision variables can be something like that:

$z(u,v) geq y_u x_v$ , $quad x,y,z in {0,1}$

and to linearize the multiplication here, I’m introducing a new decision variable $r=xy$ and these constraints has been added:

$z(u,v) geq r(u,v)$

$r(u,v) leq y_u$

$r(u,v) leq x_v$

$r(u,v) geq y_u + x_v -1$

Keeping in mind that this is for undirected graph where $r(u,v)$ is different from $r(v,u)$. How can I model this in Python using Pulp and NetworkX ?